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Zusammenfassung 


Materialien von industrieller Relevanz weisen oft eine komplexe 
Mikrostruktur auf, welche einen direkten Einfluss auf das makroskopi- 
sche Materialverhalten hat. Um diese Mikrostrukturinformationen 
bei Simulationen auf der Bauteilebene zu verwenden, haben sich 
Multiskalenmethoden etabliert. Insbesondere Homogenisierungsmetho- 
den werden häufig eingesetzt, da sie auf bewiesenen mathematischen 
Resultaten basieren. 

Der erste Abschnitt dieser Arbeit konzentriert sich auf die Charakteri- 
sierung komplexer Mikrostrukturen. Es wird die Anwendbarkeit von 
Minkowskitensoren, welche ihren Ursprung in der stochastischen Geo- 
metrie haben, untersucht. Hierbei wird insbesondere ein normalisierter 
Tensor als geeigneter Charakterisierer für Mikrostrukturen identifiziert. 
Als Grundlage der darauffolgenden Kapitel dient ein Homogenisie- 
rungsresultat für das Mumford-Shah Funktional mit periodischen 
Parametern, welches für Funktionen mit beschränkter Variation definiert 
ist. Diese Funktionen bilden auch die mathematische Grundlage für das 
Sprödbruch Modell von Francfort und Marigo. Das Homogenisierungs- 
resultat ist für den Fall einer Scherung senkrecht zur Ebene sowie unter 
Vernachlässigung der Irreversibilität auf das Francfort-Marigo Modell 
anwendbar. Erweiterungen des Homogenisierungsresultats lösen diese 
beiden Beschränkungen individuell auf und übertragen das Resultat auf 
stochastische Mikrostrukturen. Das resultierende homogenisierte Modell 
weist eine ähnliche Struktur zum ursprünglichen Francfort-Marigo 
Modell auf, allerdings mit potentiell anisotroper effektiver Steifigkeit 
und richtungsabhängiger effektiver Rissenergie. Zudem beinhaltet das 
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Homogenisierungsresultat spezifische Zellformeln fiir die effektiven 
Materialparameter, welche voneinander entkoppelt sind. 

Der zweite Teil dieser Arbeit widmet sich der Entwicklung numerischer 
Verfahren zur Berechnung der effektiven Rissenergie heterogener 
Materialien. Die Zellformel zur Berechnung dieser lasst sich als konvexes, 
jedoch nicht streng konvexes Optimierungsproblem aufstellen, welches 
den (periodischen) minimalen Schnitt durch eine Mikrostruktur be- 
stimmt. Zur Lösung dieses Problems werden eine neue Diskretisierung 
sowie FFT-basierte Lösungsverfahren mit adaptiver Parameterwahl 
entwickelt. Mithilfe dieser Strategie können Materialien von industrieller 
Komplexität untersucht werden. 

Im darauffolgenden Abschnitt wird der Fokus auf lokal anisotrope 
Mikrostrukturen gelegt. Diese erfordern eine Umformulierung des 
Zellproblems in eine anisotrope Form sowie eine Modifizierung des 
Lösungsverfahrens. Mithilfe dieser Erweiterungen können polykris- 
talline Materialien sowie Faserstrukturen mit anisotropen Fasern 
untersucht werden. 

Der letzte Abschnitt ist Untersuchungen zu verschiedenen Randbedin- 
gungen gewidmet. Die ursprünglich hergeleitete Zellformel verwendet 
Dirichlet-Randbedingungen. Untersuchungen zur Homogenisierung 
von Elastizitäts- oder Wärmeleitungsproblemen zeigen jedoch, dass 
die Verwendung periodischer Randbedingungen zu einem geringeren 
Fehler der approximierten effektiven Eigenschaften führt. Basierend 
auf diesen Ergebnissen wird untersucht, ob diese Verbesserung auf die 
Berechnung der effektiven Rissenergie übertragbar ist. 
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Summary 


Materials used in an industrial context often exhibit a complex mi- 
crostructure which directly influences the macroscopic material be- 
havior. For simulations on the component scale, multi-scale meth- 
ods may exploit the microstructural information. In particular ho- 
mogenization methods are often used due to their well formulated 
mathematical background. 

In a first step we focus on the characterization of complex microstruc- 
tures. We investigate the applicability of Minkowski tensors, which 
originate from stochastic geometry, to characterize microstructures. We 
identify in particular a normalized tensor, the quadratic normal tensor, 
as a suitable characterizer. 

The foundation of the subsequent chapters is laid by a periodic homoge- 
nization result for free discontinuity problems. These free discontinuity 
problems form the mathematical basis of the variational approach to 
fracture, known in the fracture community as the Francfort-Marigo 
model of brittle fracture. The homogenization result is directly applicable 
to a single loading step of the Franfort-Marigo model in the special case 
of anti-plane shear and without an irreversibility constraint. Further- 
more, recent extensions lift these restrictions individually and provide 
extensions to random ergodic media, indicating that the homogenization 
result holds for the general Francfort-Marigo model and non-periodic 
microstructures. The homogenized model has a similar structure as 
the original Francfort-Marigo model but with homogeneous, possibly 
anistotropic material parameters, namely an effective stiffness and an 
effective crack energy. The homogenization result includes specific cell 
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formulas for both effective properties which decouple upon homoge- 
nization. The second part of this thesis is devoted to establish numerical 
tools to compute the effective crack energy of heterogeneous materials. 
The cell formula to compute the effective crack energy may be for- 
mulated as a convex, but not strictly convex, optimization problem 
to compute the (periodic) minimum cut through a given microstructure. 
To solve this problem we investigate novel discretization methods and 
FFT-based solvers with an adaptive parameter choice. With this strategy 
we investigate microstructures of industrial complexity. 

Subsequently, we focus on locally anisotropic microstructures. This 
requires an anisotropic formulation of the governing cell formula and 
a modification of the solver. These alterations at hand, we investigate 
polycrystalline materials and fibrous structures with a distinct anisotropy 
of the fibers. 

The last section is devoted to investigating the influence of different 
boundary conditions. The original cell formula for the effective crack en- 
ergy is formulated using Dirichlet boundary conditions. Investigations 
on the homogenization of linear elasticity or heat conductivity show an 
advantage of periodic boundary conditions over the Dirichlet kind in 
the convergence rate of the apparent properties. Based on these results 
we investigate if this advantage translates to the computation of the 
effective crack energy. 
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Chapter 1 


Introduction 


1.1 Motivation and objectives 


With fracture, we denote a separation of a previously coherent material 
due to breakage of bonds between the molecules of the material. 
This breakage may occur due to crack initiation, the propagation 
of a pre-existing crack through the material or a complete failure of 
the structure. Fracture mechanics describes these different cracking 
behaviors of materials from various perspectives. The field of research 
may be broadly categorized into two principal cases, namely brittle 
fracture and ductile fracture. Brittle fracture describes a sudden crack 
propagation with almost no energy dissipation before the fracture. For 
ductile fracture on the other hand, a noticeable amount of energy has 
been dissipated before the crack propagates in the form of damage and 
plasticity. Furthermore, one distinguishes the dynamic case and the 
quasi-static case. 

The actual mechanisms of fracture and crack propagation take place over 
a broad range of length scales. On the atomic scale, the separation of a 
previously intact material is described by breakage of inter-atomic bonds. 
On the component scale, a continuum mechanical perspective is typically 
taken. However, the classical continuum mechanical toolkit faces 
difficulties in the presence of cracks, which make in the displacement 
field discontinuous. Additional difficulties arise in the presence of 
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propagating cracks where the newly formed crack surface provides an 
additional unknown. 

One particular question of interest in the field of brittle fracture 
is when cracks propagate, which is the basic motivation for failure 
criteria. Failure criteria play a crucial role in the engineering design of 
components to ensure their safety and reliability. Typical failure criteria 
involve either a stress-based aproach, are energy based or focus on stress 
intensity factors. To equip failure criteria with necessary parameters, 
crack-related material properties are required. Of particular interest in 
brittle fracture are the critical energy release rate, the fracture toughness, 
as well as the ultimate tensile strength. 

The theoretical strength of materials, derived from the strength of 
the atomic bonds in a hypothetically perfect arrangement, typically 
overestimates the actual strength measured in real-world experiments. 
This deviation is caused by a non-perfect arrangement of atoms on a 
small length scale mainly in the form of micro cracks and defects. 

In addition to flaws and micro cracks, industrial materials often exhibit a 
random microstructure due to their manufacturing process. From 
fiber reinforced composites to polycrystalline materials or porous 
media, several intended or unintended production factors influence the 
microstructure. This makes a fully experimental characterization rather 
expensive. An effective tool to incorporate microstructure information 
into material models is given by multi-scale simulation methods and, 
in particular, homogenization methods. Assuming a distinct scale 
separation, the aim of homogenization methods is to treat the material 
on the component scale as homogeneous, where effective mechanical 
quantities and material coefficients arise from a separate field problem 
on a microstructure cell. 

This thesis is devoted to the characterization of microstructures and 
the investigation of a homogenization result for variational brittle 
fracture (Francfort and Marigo, 1998). For the characterization of 
microstructures we rely on an approach originated in stochastic geom- 
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etry, the Minkowski tensors. In order to characterize microstructures 
of different sizes we establish the quadratic normal tensor as our 
primary characterizer. 

For the homogenization of brittle fracture we pursue an approach 
based on a periodic homogenization result (Braides et al., 1996) for 
the Mumford-Shah functional (Mumford and Shah, 1989), as well as 
several extensions including irreversibility (Giacomini and Ponsiglione, 
2006), random, ergodic media (Cagnetti et al., 2019) and the full 
Francfort-Marigo functional (Friedrich et al., 2022). Based on these 
theoretical works we consider a cell formula defining the effective crack 
energy. We aim to provide numerical tools to compute the effective crack 
energy building on an initial work of Schneider (2020). 


1.2 Outline and originality 


Chapter 2 is devoted to the concepts which form the basis of this work. 
We first give an overview on brittle fracture mechanics. We discuss the 
variational approach to fracture and the typical numerical treatment of 
the latter via phase-field fracture models. A second fundamental concept 
is given by microstructure characterization which is required to distin- 
guish different microstructues using simple quantities. Based on these 
characteristics, homogenization methods may be applied. Thirdly, we 
discuss phase-field fracture on complex microstructures. The governing 
equations require careful treatment and we propose an implicit solver 
strategy. Lastly, we discuss a homogenization result (Braides et al., 1996) 
for the Mumford-Shah functional (Mumford and Shah, 1989) which 
shares various properties with the Francfort-Marigo model of brittle frac- 
ture (Francfort and Marigo, 1998). We discuss the requirements for this 
homogenization result to be applicable to the Francfort-Marigo model. 
We define the effective crack energy based on a cell formula derived by 
Braides et al. (1996). Extending Schneider (2020) we identify this formula 
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as the minimum cut problem and derive an equivalent form based on 
the minimum cut/maximum flow duality derived by Strang (1983). 

In chapter 3 we focus on the characterization of digital microstructures. 
We propose the use of Minkowski tensors, which originate in stochastic 
geometry, as a tool to characterize microstructures. These tensors are 
applicable to a wide range of different shapes. In particular, we introduce 
the quadratic normal tensor as anormalized Minkowski tensor which 
provides a measure of the anisotropy of a microstructure. We introduce 
anumerical method to compute the quadratic normal tensor on large 
digital microstructures. We compare the accuracy of our approach 
for fiber reinforced composites with the usual tool to determine the 
fiber orientation tensor based on a structure tensor approach. A final 
investigation on sand grain structures demonstrates the wide range of 
applicability of the quadratic normal tensor. 

Chapter 4 is concerned with efficient algorithms to compute the effective 
crack energy on heterogeneous random microstructures. In order to find 
solutions of the cell problem derived in chapter 2, suitable discretizations 
and solvers are required. We provide extensions to the approach of 
Schneider (2020) who derived FFT-based solvers and discretizations. We 
discuss an implementation of a discretization based on a combinatorially 
consistent grid, which was introduced by Couprie et al. (2011), into the 
context of FFT-based solvers. Furthermore we rely on the alternating 
direction method of multipliers to solve the discretized problem using 
an adaptive parameter strategy. 

In chapter 5 we use both solvers and discretizations provided in the 
previous chapter and propose an extension to the anisotropic case. We 
discuss the anisotropic minimum cut/maximum flow problem and 
apply it to study anisotropic materials. 

Chapter 6 is devoted to explore the influence of the boundary conditions 
when computing the effective crack energy on random microstructure 
cells of finite size. For this purpose we use fast marching methods to 
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compute the effective crack energy in two dimensional structures which 
allows us to choose the boundary conditions freely. 

Finally, we summarize the novelties of our approach and give a final 
conclusion in chapter 7 . 


Chapter 2 


Fundamental concepts 


2.1 Basic concepts of brittle fracture 
mechanics’ 


2.1.1 Linear elastic fracture mechanics 


Modern fracture mechanics (Gross and Seelig, 2017) originated from 
the pioneering work of Griffith (1921), who postulated an energetic 
criterion for the quasi-static growth of a pre-existing crack in a brittle, 
isotropic, and elastic solid. He considered the change of the potential 
energy II with the crack area A and postulated that a crack can only 
grow whenever the energy release rate —dII/dA reaches a critical value 
y. Put differently, the crack grows whenever it is energetically more 
favorable to increase the surface energy of the crack than to increase the 
elastic energy stored in the body. 

For the case of three-dimensional isotropic elasticity and a semi-infinite 
planar pre-existing crack in an infinite medium, Irwin (1957) investigated 
the stress state at the crack tip, which has an r~ 2-singularity where 
r denotes the Euclidean distance to the crack tip. He distinguished 
three different modes which quantify the degree of singularity of the 


1 This section is based on the introductory sections of Ernesti and Schneider (2021; 2022) 
and Ernesti et al. (2023) and provides an extension thereof. 
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(a) Mode I (b) Mode II (c) Mode II 


Figure 2.1: Schematic illustration of mode I, mode II and mode IN loading. 


normal, the in-plane and the out-of-plane shear stresses, respectively, see 
Fig. 2.1 for a graphic illustration. He associated a stress-intensity factor 
to each mode and postulated that the pre-existing crack has to grow 
provided a combination of the stress intensity factors reaches a critical 
value, the so-called fracture toughness. Expressing the energy-release 
rate in terms of these stress-intensity factors, Irwin converted Griffith’s 
criterion into an investigation of the associated stress concentration at 
the crack tip (Irwin, 1962). Thus, in the case of linear elastic, isotropic 
and homogeneous brittle materials both Griffith’s energy criterion and 
Irwin’s notion of stress intensity factors coincide. In the literature, 
fracture toughness and critical energy release rate or crack resistance are often 
used interchangeably (Gross and Seelig, 2017). This work is focused on 
Griffith’s point of view and we call the energy related quantity the crack 
resistance, avoiding the term fracture toughness altogether. 

Independently, Cherepanov (1967) and Rice (1968) proposed a way to 
compute the local energy-release rate at a crack tip in terms of a contour 
integral around the crack tip, the so-called J-integral. This J-integral 
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offers an elegant way to deal with the stress singularity at the crack tip 
numerically, as it is independent of the chosen path of integration. 
Following the first concept of Irwin (1957), the stress is infinitely large 
at the crack tip. This motivated a refinement of the model to account 
for small plastic effects close to the crack tip since atomic bonds cannot 
withstand an infinite load. These small-scale plastic effects are even 
considered for macroscopic brittle material behavior. The small plastic 
region around the crack is called the plastic zone (Irwin, 1968). 

All three aforementioned pillars - Griffith’s energy release rate, Irwin’s 
stress intensity factors and techniques based on the J-integral - may be 
extended to (homogeneous) anisotropic brittle solids (Sih et al., 1965; 
Wu, 1967; Saouma et al., 1987; Williams, 1989). In the anisotropic case, 
however, some care has to be taken, as the fracture modes are not directly 
correlated to the normal and tangential jumps of the displacement 
field (Laubie, 2013). 

Linear elastic fracture mechanics was also extended to account for elasto- 
plastic effects. Dugdale (1960) investigated mode I crack propagation, 
see Fig. 2.1a, with a perfectly plastic material behavior within the plastic 
zone. In his model, the plastic zone has an elongated shape in the 
anticipated crack direction. A generalization of this model is found in 
the form of cohesive zone models (Barenblatt, 1962; Elices et al., 2002), 
which express the correlation of the crack opening and the governing 
stresses explicitly via traction separation laws (Amidi and Wang, 2017). 
The concept of classical linear elastic fracture mechanics has also 
been extended to study fatigue (Rice, 1967). Additionally, leaving 
the quasi-static domain, dynamic effects like crack branching (Katzav 
et al., 2007) have been studied by accounting for dynamical stress 
intensity factors dependent on the speed of the crack. Creep and viscous 
effects have also been included (Hollstein and Kienzler, 1988). 

Linear elastic fracture mechanics is first and foremost a tool to investigate 
when a pre-existing crack propagates, which is expressed via fracture 
criteria. Of additional interest is predicting how a crack grows and to 
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determine possible crack paths and surfaces (Chambolle et al., 2009). 
The question is whether the crack propagates in the same direction 
as before, kinks at a certain angle or, e.g., splits into several crack 
branches. The most common approaches to predict these kinking 
angles exploit the principle of local symmetry (Gol’dstein and Salganik, 
1974) or follow the postulate of maximum energy-release (Hussain 
et al., 1974). In general, these two models yield different results for the 
kinking angle (Amestoy and Leblond, 1992). Of particular interest is the 
change of the crack direction at material interfaces in the context of a 
heterogeneous material. Early attempts study crack penetration at the 
interface (Cook and Erdogan, 1972; Erdogan and Biricikoglu, 1973) or 
whether the crack is deflected at the interface (Goree and Venezia, 1977; 
He and Hutchinson, 1989). 

In addition to the aforementioned analytical methods, several con- 
tributions are concerned with computational fracture mechanics, see 
Sedmak (2018) for an overview. Equipped with numerical tools, these 
may either serve to investigate when a structure fails, or how a crack 
may propagate. In particular, in the absence of analytical solutions, 
discretization methods are required. Similar to other branches of 
mechanics, the most common technique for discretizing the involved 
equations are finite elements. Early attempts to compute stress-intensity 
factors numerically have been proposed by Chan et al. (1970) and Rice 
and Tracey (1973). However, it turns out to be difficult to resolve the 
singularity at the crack tip, since the error of a finite element solution 
correlates with the regularity of the approximated fields. For this 
purpose, enriched (Moés et al., 1999) or extended (Fries and Belytschko, 
2010) finite-element discretizations were developed which account for 
the crack surface in cracked elements by adding special ansatz functions. 
An alternative approach form cohesive zone elements (Chowdhury and 
Narasimhan, 2000), which may be incorporated into the finite-element 
mesh. In these elements specific traction-separation laws (Wimmer et al., 
2009) may be implemented to model the correlation of crack opening and 
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governing stresses. These traction separation laws are readily calibrated 
by experiments (Dastjerdi et al., 2013; Amidi and Wang, 2017). However, 
caution has to be taken, as cohesive zone elements may induce a mesh 
dependency (Rimoli and Rojas, 2015). 


2.1.2 The variational approach to fracture 


Francfort and Marigo (1998) revisited Griffiths energetic fracture crite- 
rion (Griffith, 1921) and introduced a variational approach to fracture. 
Consider a given domain Q C R, d € {2,3} with sufficiently smooth 
boundary, as well as a field of stiffness tensors C : Q — L(Sym(d)) 
a field of positive crack resistances y : Q — Ryo. We denote by 
L(Sym(d)) the field of linear operators on symmetric d x d matrices and 
explicitly account for positive definite stiffnesses. The Francfort-Marigo 
model considers quasi-static crack growth. After a fixed discretization in 
(pseudo)-time, which discretizes an increasing external load, it seeks the 
displacement field u : Q — R? and the d — 1 dimensional crack surface 
S C Q by minimizing the energy functional 


FM(u, S) = : Veu(a) : C(x) : V u(x) de +f y(z)dA (2.1) 
2 Joys s 

for each time step subjected to appropriate boundary conditions. The 
quantity V°u denotes the strain field expressed via the symmetrized 
gradient operator V*. A physical constraint is given by the irreversibility 
condition that the crack surface of the current time step must contain 
the crack surface of the previous time step. The Francfort-Marigo 
functional consists of the sum of two energies, a bulk energy and a 
surface energy. The bulk energy is quadratic in the strain field V°u and 
describes the elastic deformation of the solid. The surface term is given 
by the y-weighted surface area of the crack surface S. This additive 
decomposition into two energies allows the model to describe both 
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crack initiation and crack propagation within a single model. Upon an 
increased load, a global minimum may be found by either further elastic 
deformation of the solid, or via the formation of a crack surface. Notice, 
however, the differences compared to the formulation of Griffith (1921). 
His criterion for crack propagation is concerned with local stationary 
points. This specifically includes local minima, local maxima and saddle 
points. Solutions of the Francfort-Marigo model on the other hand rely 
on a global minimization of the functional (Chambolle and Crismale, 
2019). This global minimization is, as mentioned by Francfort and 
Marigo (1998), "not dictated by any known thermodynamical argument," 
but merely a necessity for the mathematical treatment of the model. In 
contrast to the framework of linear elastic fracture mechanics based on 
Griffith (1921) and Irwin (1957), no preexisting crack is required in the 
Francfort-Marigo model. 

The Francfort-Marigo functional (2.1) shows strong similarities with 
Mumford-Shah type functionals (Mumford and Shah, 1989) used in 
image segmentation. Consider a scalar variable v : Q — R from the 
space of special functions of bounded variation (SBV) (Braides, 1998) 
with jump set S, and the symmetric, positive definite tensor field A. A 
similar functional to (2.1) from the family of Mumford-Shah functionals 
is given by 


MS(v) = TAG - A(x) - Vv(x) dx +f y(x) dA. (2.2) 
Notice three main differences compared to the Francfort-Marigo model: 
The first difference regards the bulk energy. The Mumford-Shah func- 
tional considers a quadratic function of the gradient of v, whereas the 
Francfort-Marigo model considers a quadratic function of the sym- 
metrized gradient of u. For the special case of anti-plane shear, for 
which the displacement field u = (ux, uy, uz)’ may be expressed using 
a scalar variable v = u,, which only depends on the spatial dimensions 
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x and y, the bulk energies coincide. The second difference concerns 
the integral range of the surface part. The crack surface S serves as an 
additional objective variable, whereas the surface range in the Mumford- 
Shah functional is given by the jump set of v, namely S,. The third 
difference between the models is given by the discretization in time. 
Within the Mumford-Shah functional there is no time stepping involved. 
The Francfort-Marigo model on the other hand explicitely accounts 
for changes in an external load and minimizes the Francfort-Marigo 
functional in each time step. Furthermore, the history of the previous 
steps is taken into account via the irreversibility constraint. To sum- 
marize, in the case of anti-plane shear, the Mumford-Shah functional 
describes the relaxed form of the Francfort-Marigo model in SBV where 
irreversibility is neglected. To prove the existance of weak solutions 
of the Francfort-Marigo model in the general case, i.e., including the 
symmetrized gradient instead of the gradient as well as accounting 
for irreversibility, the solution space of generalized special functions of 
bounded deformation (GSBD) (Dal Maso, 2013) is required (Chambolle 
and Crismale, 2021). In this function space the crack surface is identified 
with the jump set of u, namely Su. The existence of strong solutions 
which treat the displacement field and the crack surface as two variables 
has been shown recently (Chambolle and Crismale, 2019). Various 
approaches to compute minimizers of the functional numerically have 
been proposed (Giacomini and Ponsiglione, 2003; Pandolfi et al., 2013; 
Schmidt et al., 2009). In particular, the Francfort-Marigo model has laid 
the foundation of the field of phase-field fracture (Bourdin et al., 2000) 
which enjoys great popularity, see Wu et al. (2020) for a recent review. 


2.1.3 Phase-field fracture 
The phase-field model of brittle fracture was introduced as a regu- 


larization of the Francfort-Marigo model, motivated by the Ambrosio- 
Tortorelli approximation (Ambrosio and Tortorelli, 1990) of the Mumford- 


13 


2 Fundamental concepts 


Shah functional (Mumford and Shah, 1989). Owing to their ability to 
nucleate cracks and to produce complex crack patterns, phase-field 
fracture models were subject to a flurry of activity (Ambati et al., 2015; 
Wu et al., 2020). The phase-field model introduces a damage variable 
d : Q — [0,1] which takes values d = 0 when the material is fully 
intact and d = 1 within a crack. The range between these values is 
given by a smeared interface of phase-field width ! > 0 which provides 
a regularization of the sharp interface which represents the crack in 
the Francfort-Marigo model. After a discretization in (pseudo-)time, 
the phase-field model in its original form is given by minimizing the 


functional 


1 208 s lal? l 2 
z079) Veu:C: Veuty + z|Ival dx 


2l 
(2.3) 


PFi(u,d) = / 


Q 


in each time step for both u and d. To ensure the irreversibility of the 
phase-field model, these minimizers are sought under the constraint that 
d'+! > d’ almost everywhere, where the index i represents the time step. 
Chambolle (2004) showed that the phase-field functional I-converges 
to the Francfort-Marigo model as / — 0. Similar to the solution theory 
for the Francfort-Marigo model, this convergence result relies on global 
minimizers. 

Several methods have been established to enforce the irreversibility 
constraint, which is necessary for the physicality of the model since it 
prevents cracks from healing. Miehe et al. (2010a) introduced a formula- 
tion based on a history field which enjoys great popularity. Another ap- 
proach (Bourdin et al., 2000; Burke et al., 2010) enforces the irreversibility 
constraint by fixing d = 1 in all points x where the damage field exceeds 
a certain threshold via Dirichlet boundary conditions. Alternatives rely 
on established methods for constrained optimization problems such 
as penalty methods (Gerasimov and De Lorenzis, 2019) or augmented 
Lagrangian methods (Wheeler et al., 2014). 
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(a) Tension (b) Compression 


Figure 2.2: Schematic illustration of tension and compression under mode I loading. 


Originating from the form (2.3), several extensions have been proposed. 
The model (2.3) does not distinguish tension and compression since the 
strain field V*u enters quadratically into the functional. This is, however, 
unphysical. Consider a pre-existing crack under mode I tensile load, 
see Fig. 2.2a. The load, if sufficiently large, causes the crack faces to 
open and the crack to further propagate. Subjected to the same load in 
opposite direction, i.e., a compressive load, see Fig. 2.2b, the material 
behavior is fundamentally different since the load forces the crack faces 
to close which results in an entirely different stress state. To account 
for this distinction of tension and compression, several models have 
been proposed. Amor et al. (2009) established a splitting for isotropic 
materials based on the sign of the strain trace. A different approach was 
proposed by Miehe et al. (2010b) based on an eigenvalue decomposition 


15 


2 Fundamental concepts 


of the strain field. Chambolle et al. (2018) proved that the approach of 
Amor et al. (2009) [’-converges to the Francfort-Marigo model with the 
additional constraint that the crack faces do not interpenetrate as | — 0. 
In the same work they showed that the method which uses the splitting 
of Miehe et al. (2010b) converges to a Francfort-Marigo type model which 
prevents crack propagation due to shear loading. Recent developments 
propose extensions to anisotropic stiffnesses (van Dijk et al., 2020) and 
multi-axial loading conditions (De Lorenzis and Maurini, 2021). 

The phase-field model (2.3) is formulated within the context of brittle 
fracture and small strain elasticity. Extensions to finite deformations 
(Hesch and Weinberg, 2014; Miehe et al., 2016) and ductile fracture 
(Ambati et al., 2016; Kuhn et al., 2016) are, however, well established. 
Furthermore, applications of phase-field models to dynamic fracture 
problems were proposed (Borden et al., 2012; Mandal et al., 2020; Ren 
et al., 2019; Weinberg and Wieners, 2022) which rely on continuum wave 
equations instead of the quasi-static case. In these models, typically 
observed dynamical effects such as crack branching (Sun et al., 2021) can 
be modeled (Hofacker and Miehe, 2012). 

Phase-field models may be used to simulate crack initiation and propa- 
gation on the component scale where the material properties are often 
modeled as isotropic and homogeneous. Moreover, models which specif- 
ically account for heterogeneities on the microscale were established 
(Kuhn et al., 2015; Chen et al., 2019; Ernesti et al., 2020). 

Additionally, several extensions accounting for anisotropic crack re- 
sistances were proposed. Francfort and Marigo (1998) already pro- 
posed a way to incorporate a direction dependent crack resistance into 
their model. Focardi (2001) established a -convergence result for an 
Ambrosio-Tortorelli approximation (Ambrosio and Tortorelli, 1990) of a 
Mumford-Shah functional with a direction dependent surface term. This 
result may be applied to the weak form of the Francfort-Marigo problem 
in case of anti-plane shear, neglecting irreversibility. Implementations of 
phase-field models with an anisotropic crack resistance were proposed 
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by Clayton and Knap (2014; 2015), who introduced a geometrically 
nonlinear phase-field model and expressed the anisotropy of the crack 
resistance via a second order tensor. Na and Sun (2018) and Bryant 
and Sun (2018) proposed a method coupling crystal plasticity with an 
anisotropic phase field fracture model. Expressing a general, anisotropic 
crack resistance via a second order tensor allows for only one weak 
direction within the plane. Nguyen et al. (2017) proposed a multi phase- 
field approach to consider more complex forms of anisotropy. Other 
approaches rely on a higher order phase-field model which accounts 
for second derivatives of the damage variable d. In this case, a fourth 
order crack resistance tensor may account for a more general case (Li 
et al., 2015; Teichtmeister et al., 2017; Kakouris and Triantafyllou, 2019; 
Ma and Sun, 2020). 

Kuhn et al. (2015) investigated the use of different degradation functions 
within phase-field fracture models. Furthermore, different approxi- 
mations of the crack energy, which establish a different interpretation 
of phase-field models were proposed (Pham et al., 2011). Phase-field 
models exhibit similarities to non-local damage models (Jirasek, 1998) 
and may be treated as such, instead of as the regularization of a free- 
discontinuity problem. In this interpretation, the phase-field width 
| serves as a material parameter (Kuhn, 2013) which influences the 
ultimate tensile strength of the material based on a one-dimensional 
analysis of the model. Additionally, length-scale insensitive methods 
were proposed (Wu and Nguyen, 2018) which are referred to as cohesive 
phase-field fracture. 
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2.2 Characterization of digital microstructures’ 


2.2.1 Objectives 


Materials of industrial interest often display a heterogeneous underlying 
microstructure which influences the macroscopic material behavior. Let 
us consider an injection molded structural component made from a 
polymer material. For improving the stiffness of the component while 
keeping its lightweight potential, filler materials, such as glass fibers, 
are commonly added to the polymer melt before being modeled. As a 
result, the component's material is given by a fiber reinforced polymer 
which has a complex microstructure determined by the arrangement 
of the fibers. Often, both polymer and glass are considered isotropic. 
However, the macroscopic material behavior is strongly dependent on 
the fiber arrangement on the microscale which may induce a distinctly 
anisotropic material behavior on the component scale. 

Due to the infinite number of possible arrangements on the microscale, 
an exhaustive experimental characterization of a composite’s microstruc- 
ture dependent anisotropic behavior may prove challenging. Hence, 
homogenization-based multi-scale methods have been developed to 
simulate the material behavior by explicitly taking the microstructural 
information into account, see Matouš et al. (2017) for a recent overview. 
These homogenization techniques compute the effective response of the 
heterogeneous material, taking the material behavior of the constituents 
and the microstructure into account. Therefore, the microstructure has 
to be quantified in terms of suitable data, which is where microstructure 
characterization comes into play. 


? This section is based on the introduction of Ernesti et al. (2022). Changes to the text have 
been made in order to include it into the structure of this work. 
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2.2.2 State of the art 


A common image-based microstructure-characterization method is scan- 
ning a material sample via micro-computed tomography (u-CT) (Schla- 
ditz et al., 2017; Cnudde et al., 2009). After tomographic reconstruction, 
the local mass density of a material is determined and stored as 3D 

voxel data. In case of a two-phase material and after some processing, 
this voxel data may be interpreted as the characteristic function of the 

microstructure, i.e., the function which attains the value 1 for one phase, 
and the value 0 for the complementary phase. Correctly segmenting 

u-CT scans requires a certain contrast in the absorption rates of the 

constituents to be applicable, for instance for porous media or for a 

variety of composite materials. 

The mechanical behavior of composite and porous materials is strongly 

influenced by the volume fractions of the phases. If the characteristic 

function is accurately resolved by the „CT-scan, the volume fraction 

can be computed accurately by numerical integration. With the volume 

fraction at hand, bounds that predict the possible range of effective 

elastic and thermal material properties may be established, see Voigt 

(1889) and Reuss (1929). However, for a high material contrast, these 

bounds span a wide range and hence provide limited information. 

For higher accuracy, additional information is required, see Torquato 

(2002) for an overview. For instance, n-point correlation functions (Brown, 
1955; Torquato and Stell, 1982) provide suitable additional information. 
Their applicability for anisotropic materials, however, is limited due to 

the high associated computational effort, see Eriksen et al. (2004). 

For the class of fiber-reinforced composites, specific microstructure- 
characterization techniques have been established. In addition to the 

fiber volume-fraction, common characteristics include the fiber aspect- 
ratio and fiber-orientation tensors of second and fourth order (Kanatani, 
1984; Advani and Tucker, 1987), see for instance Müller and Böhlke 

(2016). A variety of methods for computing fiber-orientation tensors 
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based on volumetric images has been established (Robb et al., 2007; 
Daniels et al., 2007). A common approach is based on the so-called 
structure tensor (Krause et al., 2010). 

For porous structures, different microstructure characteristics are of 
interest. For instance, the tortuosity (Neumann et al., 2019) and 
chord-length distribution (Matheron, 1975; Torquato, 2002), as well 
as the pore-size distribution (Kate and Gokhale, 2006) are investigated. 
These measures are primarily responsible for the effective (isotropic) 
permeability of the porous medium in question. 

Polycrystalline materials require a different approach. Typically, the 
grains differ only in their crystalline orientation, but have identical 
absorption rates. Hence, u-CT scans are of limited use. Instead, 
for reconstructing the 3D-microstructure of polycrystalline materials, 
focused ion beam - scanning electron microscopy (FIB-SEM) (Bansal 
et al., 2006; Groeber et al., 2006; Zaefferer et al., 2008) or electron 
back-scattering diffraction (EBSD) (Korte et al., 2011; Adams and 
Olson, 1998; Larsen et al., 2002) are preferred. Primary microstructure 
characteristics of polycrystalline materials are the grain-size distribution 
(morphological texture) (Döbrich et al., 2004) and the orientation 
distribution (crystallographic texture) (Bunge, 1982; Böhlke, 2006; Böhlke 
and Lobos, 2014; Junk et al., 2012; Bohlke et al., 2010). 

From a theoretical point of view, most materials undergoing a manufac- 
turing process are influenced by stochastic factors, for instance due to 
slight variations in the composition or the seemingly chaotic behavior 
of the processing condition as a result of a high sensitivity to initial 
and boundary conditions. Still, the experimentally determined effective 
properties of such composites are often surprisingly deterministic. 
These observations may be formalized by the theory of stochastic 
homogenization (Kozlov, 1979; Papanicolaou and Varadhan, 1981) which 
forms the basis for modern computational multi-scale methods. We refer 
to section 2.4.2 where we discuss these homogenization methods in the 
context of brittle fracture. 
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2.3 Phase-field fracture on heterogeneous mi- 
crostructures’ 


In this section, we investigate the classical phase-field fracture problem 
on heterogeneous microstructures. Classical here means that the phase- 
field functional we consider uses both a quadratic damage degradation 
function and a quadratic damage penalty term. We derive the governing 
equations in section 2.3.1. We introduce several energy splittings from 
the literature (Amor et al., 2009; Miehe et al., 2010b). Additionally, we 
follow Boeff et al. (2015) to adjust the phase-field model for heteroge- 
neous material properties. 

To solve the governing equations we rely on FFT-based solvers and 
discretizations on regular voxel grids, as well as a fully implicit scheme 
which is first order in time. To be more precise, we use a fast gradient- 
method for the staggered treatment of the fully-coupled nonlinear phase- 
field fracture system. Fast gradient-methods add an inertial term to 
classical gradient schemes, dramatically accelerating the latter, and 
exhibit good performance also for non-convex problems of interest, 
with recent research triggered by machine learning applications. The 
theoretical reasons are still under investigation (Jin et al., 2018), but it 
is commonly believed that the inertial term helps navigating through 
non-convex energy landscapes more rapidly, ignoring local stationary 
points and flat ravines. Here, we rely on the heavy-ball method, first 
introduced by Polyak (1964). 

We investigate phase-field fracture on the complex microstructures of 
fiber reinforced composites in sections 2.3.3 and 2.3.4. 


3 This section is based on Ernesti et al. (2020) which introduces a novel solver for phase- 
field fracture on heterogeneous microstructures. This solver was first presented in 
my masters thesis (Ernesti, 2018) and the publication (Ernesti et al., 2020) provides an 
extension thereof. We exclude the detailed discussion on the solver and its performance 
and provide only a selected choice of numerical examples which serve our purpose here. 
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2.3.1 Continuous model under investigation 


Let Y = (0, L1] x [0, L2] x [0, L3] be a rectangular unit cell. Suppose two 
heterogeneous elastic strain energy potentials Y+, Y~ : Y x Sym(3) > R 
are given, where Sym(3) denotes the space of symmetric two-tensors in 
R. ~~ measures the stored elastic energy per unit volume which is not 
degraded, whereas ı)* is responsible for the part which degrades as a 
result of damage. These potentials shall sum to the total (local) stored 
elastic energy density y% : Y x Sym(3) — R, i.e., the equation 


p(z, €) = wt (a,e) + (2, €) (2.4) 


shall hold for all x € Y and £ € Sym(3). For the present manuscript, 
we set d(x,e) = ġe : C(x) : e with a heterogeneous linear elastic 
tensor field C. 

Let furthermore a subvolume Y C Y be given, containing all microscopic 
points where the material can be damaged. (In particular, Y may be 
equal to Y.) Denote by x the characteristic function of Y, i.e., the function 
which equals unity on Y and vanishes elsewhere. 

For a given macroscopic strain E € Sym(3), the functional under consid- 
eration reads 


2 
F(u,d) = f f(x, d(£)jyt(£,£) +h (x, €) + Ge (5 + ; Iva?) dx, 
i (2.5) 
a function of the periodic displacement fluctuation field u: Y + R? and 
the periodic damage field d : Y — R. Here, e = E + V*u denotes the 
total strain involving the symmetrized gradient operator V°, 


F(a, d) = ko + (1 — ko)(1 — x(x)d)? (2.6) 


is the damage degradation function involving a relative residual stiffness 
ko> 0 and G. > 0 stands for the critical energy release rate according 
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to Griffith (1921). The latter is assumed to be homogeneous for now as 
we specify the material points that may damage via the characteristic 
function x. The length parameter ! > 0 describes the width of the 
transition zone between broken and intact material, and is decisive for 
crack nucleation. Before continuing, we discuss the types of splittings 
(2.4) we will consider in this section, where we assume C to be isotropic, 
i.e., its application may be written as 


C(x) : € = 2u(z)e + A(x) tr(e) Id = 2u(x)e' + K(x)tr(e) Id, 


(2.7) 
e€Sym(3), zeY, 


for positive functions A, u, K : Y + R, which encode the heterogeneous 
Lamé’s constants and the compression modulus, and the deviatoric 
strain e’ = e — %tr(e) Id. In this isotropic case we distinguish three 


different energy splittings: 


1. Trivial splitting: no distinction between tension and compression 


wt (x, e) = Le :C(x):e, Y~ (x,€)=0, xeyY, ee€Sym(3). 
(2.8) 
2. The spherical splitting of Amor et al. (2009) 


wt (z,e) = ulz)lle’||? + 2A (x) (tr(e))4, 
w(x, €) = 2X(x)(tr(e))?, (2.9) 
zeY, e€ Sym(3), 


where (-)+ denote the McCauley brackets 


(-)4 =max(0,-) and (-)_ = min(0,-). 
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3. The eigenvalue splitting of Miehe et al. (2010b) 


U" (x, e) = u(x) ({e1)4 + (€2)% + (ed$) + A(a)(tr(e)) i, (2.10) 


zeY, eeSymß), 


<A 


where c; (i = 1, 2, 3) denote the eigenvalues of the symmetric tensor £. 
The last two splittings enable the model to distinguish tensile and 
compressive loading. More precisely, Chambolle et al. (2018) proved 
that the phase-field model using Amor’s splitting T-converges for | — 0 
to the Francfort-Marigo model with a constraint of non-interpenetration. 
In the same paper it was shown that the splitting of Miehe et al. (2010b) 
converges to a Francfort-Marigo model which restricts crack opening to 
shear loading, which is typical, for instance, for concrete. 
Suppose that js and A are essentially bounded from above. Then, for 
all three splittings considered, the phase-field functional (2.5) is well- 
defined on the product Banach space V x Z, where V contains all periodic 
H' displacement vector fields with mean zero and Z is the intersection 
of periodic scalar-valued H+ damage fields with L(Y). The restriction 
to essentially bounded damage fields ensures that the first term in (2.5) 
is integrable. 
It can be shown, see Burke et al. (2013), that F is Gäteaux-differentiable 
on V x Z with directional derivative 


DF(u, d)|v, z| = D,F(u,d)[v] + DaF (u, d)[z] (2.11) 
for u,v € V,d,z € Z. The directional derivatives are given by 


+ = 
toe + er :Vivde, e=E+ Vu, 


(2.12) 


D,F (u, d)[v] = / 


Y 
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where we, for simplicity of notation, suppress the x-dependence, and 


DaF (u, d)[z] = F 


PrN jyt (e)z + Ge E +1Vd- v,| dz. (2.13) 
We say that (u, d) is a critical point of F if both D,,F(u,d) = 0 in V’ and 
DaF(u,d) = 0 in Z’ hold. Integrating by parts, it is readily seen that 
such a critical point (u, d) satisfies both the balance of linear momentum 


div (o) (c,d) = 0 (2.14) 
with the stress 
a(e,d) = f(d) ad (e) + eo), e=E+V’u, (2.15) 


and an advection-diffusion-type equation for the damage variable 


o= laute) + Ge ($ -radl 


which can, for later reference, be equivalently rewritten in the form 
(l+a)d—PAd=a (2.16) 


involving the function a : Y — R determined by the equation 


ala) = 2 (1 ko)x(0) Y+ (2, e(a). (2.17) 


To conclude the discussion of the continuous model, several remarks are 

in order. 

1. The model we presented does not account for irreversibility, i.e., it 
is purely elastic. Such a model may be reasonable for monotonic 
loading, where, in addition, no local unloading occurs. To account for 
irreversibility, several strategies have been developed. Classically, the 
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constraint d > 0 is enforced, and solved by, for instance, augmented 
Lagrangian methods (Wheeler et al., 2014), primal-dual active-set 
methods (Heister et al., 2015) or penalty formulations (Gerasimov 
and De Lorenzis, 2019). 

An alternative approach is based on the usage of a suitable history- 
field, as pioneered by Miehe et al. (2010a). More precisely, for a given 
(pseudo-)time-dependent strain loading 


E: [to, ty] = Sym(3), 


the history-field based model seeks the displacement field and the 
damage field 


u : (to, ti] x Y —> R? and d: [to, t1] x Y > R, 
s.t. the balance of linear momentum (2.14) is satisfied, for 
eft, x) = E(t) + V°u(t, x), 


together with the advection-diffusion-type equation (2.16) for the 
damage variable with the forcing term 

a(t, x) = a —ko)x(x) sup wWt(z,e(r7,2)). (2.18) 
Ge to<t<t 

Thus, the entire history of the elastic energy density is taken into 
account, in contrast to equation (2.17). From an implementation 
point of view, using either (2.17) or (2.18) does not make much of a 
difference, because equation (2.16) needs to be solved in both cases. 
However, the resulting models are different, see (Ernesti et al., 2020, 
section 4.3.6). 


. The model we presented is quasi-static, and allows for brutal fracture. 


However, some care has to be taken, as brutal fracture involves cracks 
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propagating with infinite speed. This violates physical experience, 
where cracks propagate with finite speed. To account for finite-speed 
crack-propagation, Ginzburg-Landau-type equations (Miehe et al., 
2010b; Kuhn, 2013) replacing (2.16) may be used. More precisely, 
a multiple of d is added to the equation (2.16). Upon a backwards 
Euler time-discretization, the resulting equation is of similar type 
as equation (2.16) and may be treated by the same techniques, see 
Chen et al. (2019), section 2. However, for the sake of brevity, we will 
consider the unregularized form (2.16). Investigating this choice 
has the added benefit that our computational techniques do not 
deteriorate for vanishing viscous regularization (or, equivalently, 
infinite phase-field mobility parameter). 

. The reader may wonder about the particular form of the energy 
functional (2.5) that we consider. In particular, the appearance of 
the set Y C Y (and its associated characteristic function) are non- 
standard. 

Most of the models and computational strategies for phase-field 
fracture have been developed with homogeneous materials in mind. 
Consider, for the sake of exposition and / > 0 , the functional 


a of 
Gud) = f ale dyt (ae) +e) +9] + 5 VAP] de 
Y 
(2.19) 
with g(a, d) = ko + (1 — ko)(1 — d)? and a heterogeneous crack resis- 
tance y : Y — (0,00). The corresponding Euler-Lagrange equations 
consist of the balance of linear momentum (2.14) (with g replacing f 
in (2.15)) and the phase-field equation 
o 
0= ne d)yt (e) + — l div (yVd). (2.20) 
Thus, in contrast to the phase-field equation used in this work (2.16), 
equation (2.20) involves gradients of the crack resistance y, compli- 
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cating the numerical treatment. 

Instead, we chose to work with a homogeneous crack resistance Ge, 
but restrict the domain, where cracks may occur, to the subset Fer 
Thus, we consider Y asa homogeneous elastic-brittle material with 
internal variable field d, whose evolution is governed by the non-local 
PDE (2.17). 

Clearly, our framework can be extended to heterogeneous crack 
resistances as follows. Suppose y is piece-wise constant, i.e., might be 
written as y = >, G.,kXr, Where the G., are positive and the sets 
corresponding to x; form a non-overlapping partition of Y. Then, 
the natural extension of (2.5) is the functional 


Fand... de) = [ fla,dı,:.-‚de)b*(x,e) + Y7 (ze) 
Y 
x d | 
[8 + 2 va? 
+9 Gen [FF + y Vael] ae 


involving K distinct periodic damage variables dı,...,dx, and 
fE reads 


2 
f* (a, di,...,dx) = ko + (1— ko) )(1- orators a) 


Thus, also K different phase-field equations need to be solved for 
this approach. For this section, we restrict to a single brittle phase 
Y, characterized by a single crack resistance Ge, and consider the 
complement of Y in Y to deform elastically. 

The principal advantage of our approach is that, for the phase-field 
equation (2.16), the coefficients in front of the derivatives are homo- 
geneous, implicating a simplified treatment by FFT-based methods 
compared to the formula (2.20). 


2.3 Phase-field fracture on heterogeneous microstructures 


4. Due to the x-term in (2.17) and (2.18), a is concentrated on the chosen 
region Y. Similarly, f(x,d) = 1 for x ¢ Ý. Thus, d has no direct effect 
on the mechanical behavior of the composite outside of Y. 

Formally setting l = 0 in (2.16) shows that d vanishes outside of Y. For 
positive l, the damage field may also be positive in Y\Y. However, 
for small l, d will be small in Y\ř, as well. Rather, this phenomenon 
is closely tied to computing the energy dissipated by creating a crack, 


ae 1 
— +- ||Vdl|?| d 
fo [+ giver] «. 


as a part of the functional (2.5), correctly. We have learned this ap- 
proach from Boeff et al. (2015), who call it "non-local phase damage". 


2.3.2 Numerical setup 


In order to solve the governing equations (2.14) and (2.16), suitable 
discretizations and solvers are required. In the work at hand we rely on 
an FFT-based solver framework and discretization methods on regular 
voxel grids. We integrated this into an in-house FFT-based computa- 
tional micromechanics solver, based on Python 3.7 with Cython exten- 
sions. This in-house code is strain-based, i.e., the primary variable is the 
(compatible) strain field instead of the displacement field. Furthermore, 
we consider periodic boundary conditions. Thus, we seek a periodic 
strain field with prescribed mean value E. We rely on a staggered grid 
discretization (Schneider et al., 2016) to discretize the strain field. For 
the damage variable we rely on the original Moulinec-Suquet discretiza- 
tion (Moulinec and Suquet, 1994; 1998). 

In each load step we solve the balance of linear momentum using an ac- 
celerated gradient descent method, namely the heavy-ball method, first 
introduced by Polyak (1964). In each iteration we solve the advection- 
diffusion type equation for the phase-field variable using gradient de- 
scent. I introduced this approach in my master thesis Ernesti (2018) and 
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for a detailed discussion I refer to (Ernesti et al., 2020, section 3). 

To change between load steps, affine extrapolation is used for the strain, 
as described by Moulinec and Suquet (1998), and the last converged 
damage field is taken as initial guess for the succeeding step. 

For the mechanical sub-problem, we use the convergence criterion pro- 
vided in (Ernesti et al., 2020, equation (3.4)), with a prescibed tolerance 
tol. In section 2.3.3 we chose tol = 107° and in section 2.3.4 we chose 
tol = 1074. The damage sub-problem is terminated upon satisfying the 
criterion in (Ernesti et al., 2020, equation(3.17)) with tol = 10~°, similar 
to Chen et al. (2019). 

The average stress X is computed by volume averaging 


1 
Dı= (ly = [rd 


of the individual components of the stress tensor. 
All computations were performed on a desktop computer with a 6-core 
Intel i7 CPU and 32GB RAM. 


E in GPa v Gein N/mm 
Polyamide 3.45 0.39 0.1 
E-glass 72.00 0.22 g 


Table 2.1: Material parameters of polyamide and E-glass (Ernesti, 2018) 


hin um | lin um ko 
2 6 1074 


Table 2.2: Numerical parameters 
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2.3.3 Phase-field fracture for a continuous-fiber rein- 
forced brittle composite 


0.0 


EEE 


SEE 


(a) #1 (b) #2 (c) #3 (d) #4 


(e) #5 (£) #6 (g) #7 (h) #8 


Figure 2.3: Damage fields and crack paths upon final failure for eight different microstruc- 
tures. (Ernesti et al., 2020) 


In this section, we investigate the phase-field fracture behavior of 
continuous-fiber reinforced composites. This allows us to rely on 2D 
structures. We consider the inclusions linear elastic and the matrix 
material brittle, see Tab. 2.1 for the material parameters. The numerical 
parameters, i.e., the pixel length h, the phase-field width ! and the 
remaining stiffness ko are listed in Tab. 2.2. These parameters have 
been selected carefully in (Ernesti et al., 2020, section 4). We consider a 
stochastic arrangement of the fillers. The fiber volume-fraction is chosen 
as 45%, with 16 fibers placed within the volume, using the Torquato-Jiao 
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algorithm (Torquato and Jiao, 2010). We rely upon images of a resolution 
of 256 x 256 pixels. 


m m m 5 m 
+# lx# 204 344 4 FF 9 #6 FT #8 


100 F 


Mee in MPa 


l l l l 
og 0.2 04 06 0.8 1 1.2 1.4 


Ex, in % 


Figure 2.4: Stress-strain diagram upon uni-axial extension, on 8 different microstructures, 
see Fig. 2.3. (Ernesti et al., 2020) 


Eight different microstructures were generated, see Fig. 2.3. These 
microstructures were subjected to uni-axial strain loading in x-direction 
and 0.1%-steps until complete failure of the structure. The resulting 
stress-strain diagrams for the stresses perpendicular to the fiber direction 
are shown in Fig. 2.4. Only small scatter is observed prior to failing. Also, 
the post-critical stresses are identical for the different microstructures. 
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structure | Etin% | Xmax in MPa 
#1 1.3 114.10 
#2 1.3 113.82 
#3 13 113.76 
#4 1.3 114.66 
#5 1.3 114.44 
#6 1.2 112.27 
#7 13 115.27 
#8 13 114.58 

uEo 1.29 + 0.03 | 114.11 # 0.83 


Table 2.3: Strain and stress at failure on the 8 microstructure realizations, see Fig. 2.3. 


Seven of the eight microstructures fail at 1.3% axial strain, see Tab. 2.3. 
Also, the maximum effective axial-stresses are very similar, at about 
114MPa and a standard deviation below 1MPa. Taking into account 
the differences in the resulting crack paths, see Fig. 2.3, this may be 
surprising. 
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2.3.4 Phase-field fracture for a short-fiber reinforced brit- 
tle composite 


/ 
y 


(a) Original microstructure (b) Crack at Ex. = 0.9% 


Figure 2.5: Microstructure and crack surface for uni-axial strain loading in x-direction. 
(Ernesti et al., 2020) 


As our final example for phase-field fracture, we consider a short-fiber 
reinforced composite with distributed fibers. The material parameters 
are identical to the previous section, listed in Tab. 2.1. 

Aligned fibers with a length of 320jm and a diameter of 324m were 
dispersed in a periodic unit cell with dimensions 256 x 256 x 2048 um?, 
up to a fiber volume fraction of 18%. Both the fiber diameter and the 
cell dimensions were chosen carefully to ensure that the numerical 
parameters in Tab. 2.2 can be used, as well. 

The resulting microstructure, see Fig. 2.5a, was generated by the Sequen- 
tial Addition and Migration algorithm (Schneider, 2017) and discretized 
by 128 x 128 x 1024 voxels. 


34 


2.3 Phase-field fracture on heterogeneous microstructures 


—— Err —>— Pyy Der 


—— Err I Yyy > Le, 


120 T T I 
100 F J 100 F 7 
80 F A 80 F J 
£ £ 
Z= 60+ | | 0) J 
A A 
40 F | 40 F Si 
20 F J 20 F 2) 
0 1 L L 0 i L L L 1 L L 
0 0.2 0.4 0.6 0.8 0 02 04 06 08 1 12 14 16 
Ezz in % Eyy in% 
(a) Stress-strain diagram for an extension in (b) Stress-strain diagram for an extension in 
x-direction y-direction 


Figure 2.6: Stress-strain diagrams for extensions in x- and y-directions. (Ernesti et al., 2020) 


The structure was subjected to uni-axial strain loading, successively 
increasing the level by 0.1% per steps. The structure failed at E,. = 0.9%, 
and the corresponding stress-strain diagrams are shown in Fig. 2.6a. The 
stress-strain diagram is linear elastic prior to failure. Also, in the diagram, 
the stresses transversal to the fiber directions are indistinguishable. 
The crack upon failure is shown in Fig.2.5b. The crack surface forms 
around the fibers, showing matrix cracking and fiber pull-out. Further- 
more, we conducted a similar experiment as before, but with loading in 
Eyy-direction. The structure failed after 16 0.1% steps. The stress-strain 
diagram is shown in Fig. 2.6b. In contrast to loading in fiber direction, a 
damaging effect is seen in the stress-strain curves. Furthermore, upon 
loading, the symmetry between Uz, and %,, is slightly broken. The 
crack upon failure is shown in Fig. 2.7b and runs entirely through the 
matrix. A transverse view of the crack is shown in Fig. 2.7d. The crack is 
straight in the x-z-plane, but avoids the individual fibers. 
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(a) Original microstructure (b) Crack at Ey, = 1.6% 


(d) Crack at Eyy = 1.6% (transverse view) 


Figure 2.7: Crack surface and stress-strain diagram for uni-axial strain loading in y- 
direction. (Ernesti et al., 2020) 
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2.3.5 Conclusions 


In this section, we presented an approach to compute phase-field frac- 
ture on heterogeneous microstructures. We discussed the governing 
equations which rely on a strategy by Boeff et al. (2015) to deal with 
heterogeneous material parameters. We solved these equations using 
FFT-based fast gradient solvers and discretizations on a regular grid. 
Finally, we investigated phase-field fracture on complex microstructures 
of fiber reinforced composites. 

The strategy presented here follows an approach often conducted in 
FFT-based micromechanics for the homogenization of complex materials, 
see Schneider (2021a) for an overview: A cell formula for the balance 
of linear momentum on a stochastic microstructure cell is conducted. 
This cell is considered large enough in order to be representative for the 
microstructure of the considered material. Using FFT-based solvers for 
this cell problem and extracting the local fields allows to compute the 
effective stress-strain relation via averaging. For linear elastic materials 
this allows to compute an effective stiffness which may then be used 
for simulations on the component scale. For hardening type materi- 
als, upscaling approaches may be conducted (Gajek et al., 2021). For 
softening type materials, however, which includes phase-field fracture, 
this relation of effective stresses and strains provides limited informa- 
tion (Gitman et al., 2007). Hence, the question remains how to extract 
effective properties in the context of fracture mechanics which result from 
a homogenization result. 
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Figure 2.8: Schematic of a crack increment in a microstructured material (Ernesti and 
Schneider, 2021, Fig. 2). 


2.4 Homogenization methods for brittle frac- 
ture’ 


2.4.1 Objectives 


As we pointed out in section 2.2, materials of industrial interest 
often exhibit a stochastic microstructure which strongly influences 
the macroscopic material behavior. To perform mechanical simulations 
on the component scale by taking the microstructure into account, 
multi-scale methods are often conducted. These multi-scale methods 
are well established for hardening material behavior, i.e., a monotone 
stress-strain relation, see Matouš et al. (2017) for an overview. In 
this context, homogenization methods (Papanicolaou and Varadhan, 
1981; Kozlov, 1978) play a key role. Building upon a well-defined 


* This section is based on the introductory sections of Ernesti and Schneider (2021; 2022) 
as well as section 1 and 2 of Ernesti et al. (2023) and provides an extension thereof. 
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mathematical foundation, these methods aim to establish effective 
models based on the arrangement of the microstructure and the 
governing material behavior of the different phases on the microscale. 
Computational approaches typically rely on simulations on finite cells 
which represent the underlying microstructure. If these are chosen 
sufficiently large, these are called representative volume elements 
(RVEs) (Drugan and Willis, 1996; Kanit et al., 2003). 

These well established homogenization and multi-scale methods face 
difficulties when leaving the realm of hardening-type materials. In the 
presence of a propagating crack the stress-strain relation is no longer 
monotonic and thus softening of the material occurs. One example is 
given by Gitman et al. (2007) who performed numerical simulations 
using a non-local damage model. They showed that in the post-peak 
loading regime, the average stress response tends to zero as the size 
of the volume element goes to infinity. In particular, the strategy of 
computing effective quantities via simulations on representative volume 
elements — which works well for hardening-type materials — leads to 
impractical results in case of softening. 

Classically, homogenization relies on a distinct scale separation. More 
precisely, the macroscopic displacement or stress fields should vary 
slowly compared to the local fields on the microscale. In the classical 
linear elastic fracture mechanics setting and for an evolving crack, 
the stress singularity at the crack tip typically prohibits such a scale 
separation. Indeed, in classical linear elastic homogenization, the 
displacement field is split into a smooth macroscopic and a highly- 
oscillating microscopic part. Upon homogenization, the relationship 
between the macroscopic and the microscopipc displacement fields 
is simplified to a one-way coupling. This permits the transfer of 
information from the microscale to the macro scale via the effective 
stiffness. In the case of an evolving crack, the displacement field is 
no longer smooth at the crack tip, even for a homogeneous medium. 
In particular, the classical macro-micro decomposition, which was 
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successful for linear elastic homogenization, loses its promise. Moreover, 
the evolution of the crack needs to be accounted for at both scales. 
From another perspective, let us consider a non-local damage model. 
The non-locality is necessary to obtain mesh-independent results for a 
corresponding finite-element model. Thus, two scales are present in such 
a multi-scale non-local damage model, the typical scale of heterogeneity 
and the length scale of the nonlocality. In upscaling, the scale of 
heterogeneity is small, and we wish to pass to the limit of vanishing 
heterogeneity size. On the one hand, if we fix the non-local length 
scale, the nonlocallity would exceed the size of the heterogeneities upon 
homogenization. On the other hand, tying the non-local length scale 
to the size of the heterogeneities would recover the mesh-dependence 
of the model upon homogenization, and thus rendering the procedure 
illegitimate. Nevertheless, computations of phase-field fracture on 
microstructures may be pursued, see section 2.3. Furthermore using 
non-local damage models (Boeff et al., 2015; Berthier et al., 2014) to 
simulate crack propagation on microstructures are well established. Un- 
fortunately, connecting these results to macroscopic material properties 
appears challenging as the apparent stresses for softening materials are 
inherently size-dependent, see Gitman et al. (2007). 

The actual question is: Which model is suitable for the macroscopic 
modeling of fracture mechanics that results from the homogenization 
limit of a model on the microscale? An illustration of this procedure 
is given in Fig. 2.8. In the following section we identify the Francfort- 
Marigo model of brittle fracture (Francfort and Marigo, 1998) as a 
suitable candidate. 


2.4.2 Homogenization of the Francfort-Marigo model 


The strategy pursued in this thesis to homogenizing brittle fracture has 
been proposed by Schneider (2020) and is based on a homogenization 
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result by Braides et al. (1996) for free discontinuity problems. For a given 
domain 2 C R¢ we consider periodic functions f : Q x R? — [0,0o) and 
g:Qx Rx $41 — [0, 00) which, for constants c1, c2,¢3, ca > 0 satisfy 
the growth conditions 


ale < f(x, €) < ceo(1 + Itl?) YEER? and 


(2.21) 
c3(1 + |z|) < g(&,2,n) < ca(14+|z|) Yn € St Vz ER. 


Furthermore, we consider the family of functionals 


B,,(v) =| 5 (2.vo) d+ f g (Zne) dA. (2.22) 


Here, [v] = vt —v~ denotes the jump of v and n, denotes the unit 
normal to the jump set of v, S,. The length parameter 7 represents 
the periodicity of the functions f and g, see Fig. 2.9a for an illustration. 
We signify f as the bulk term of the energy functional (2.22) and g as 
the surface term of the functional. Braides et al. (1996) proved that this 
series of functionals converges to the homogeneous functional 


hom/(,,) — hom ) hom Ny ; 
Brom (y) IR; art] g (om)dA (2.23) 


as 7 — 0 in the sense of T-convergence. The expressions f?™ and 
g*™ denote the homogeneous bulk and surface terms. Notice that the 
heterogeneous and the homogeneous functionals are very similar in 
their structure. 

Additionally, Braides et al. (1996) provide specific cell formulas for both 
effective terms. Consider the infinite periodic continuation of the domain 
Q and a finite cuboid cell [0, L]? within this infinite domain. For the 
effective bulk function the governing cell formula, equipped with zero 
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(a) Periodic microstructure (b) Random microstructure 


Figure 2.9: Schematic illustration of a periodic microstructure with period n and a 
stochastic microstructure with characteristic length n. 


Dirichlet boundary conditions on the boundary of the cube, reads 


PaE = ming, A pp EVT ew 
This is a classic formula for the homogenization of bulk energies. Sup- 
pose, for instance, f(z,£) = ¿T A(x)é for some symmetric, positive 
definite matrix field A. In this case, the cell formula (2.24) may be 
interpreted as the homogenization of the static heat equation with A 
denoting a heterogeneous heat conductivity tensor and € denoting the 
gradient of the temperature field. The resulting effective bulk term 
provides a formula for the effective heat conductivity A“ (Jikov et al., 
1994; Bakhvalov and Panasenko, 1989). Alternatively to the Dirichlet 
boundary conditions used in (2.24), periodic boundary conditions are 
also suitable since the effective bulk energy does not depend on the 
boundary conditions (Sab, 1992). Furthermore, using periodic boundary 
conditions for a periodic homogenization problem permits to rely on a 
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yı 


O D 


Figure 2.10: Visualization of the computation of the effective crack energy. The cube LQn 
is placed into the periodic structure. Depending on the material contrast y2/y1 either the 
green or the red path is favored. The image is a slight modification of Fig. 1 in Ernesti et al. 
(2023). 


single period of the structure. For this case, well established tools for 
solving the periodic cell problem are available (Moulinec and Suquet, 
1994; 1998; Dorn and Schneider, 2019). 

For the effective surface term, Braides et al. (1996) proceed as follows. 
Consider a rotated cuboid LQ,, with its eı axis aligned with the normal 
n and edge length L placed in the infinite periodic continuation of our 
domain Q, see Fig. 2.10. The formula for the effective surface term reads 


hom( yn) = lim inf 


g Loowvev L 


1 
1 / g(x, [v], mv) dA (2.25) 
LQnOSv 
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with the solution space 


Y= fu € SBV(LQ,,), Vv = 0a. e., 
(2.26) 


zxz- n>0 


Vz ondLQn \ . 


0O,2-n <0 


and z € R quantifying the jump of v at the boundary. This minimization 
problem is fundamentally different to the problem associated with 
the bulk term. In (2.25) we seek a function in the space of SBV with 
minimal g-weighted integral over its jump set. In other words, we 
seek a g-weighted minimum surface cutting through the rotated cube 
LQ, with mean normal n fulfilling the boundary conditions specified 
in the set V. Thus, we refer to (2.25) as the minimum cut problem. A 
visualization of this procedure for a periodic structure of inclusions 
is shown in Fig. 2.10, where we set g(-,z,-) = y2 in the inclusion and 
g(-, £,+) = Yı in the surrounding area. In this example we distinguish 
two cases, depending on the material contrast 72/7. For a large contrast, 
the minimum surface will avoid the inclusion, illustrated by the red 
curve, for a small contrast it will cut the inclusion in a straight line. 
Similar to the bulk term, formula (2.25) is expressed via an infinite 
volume limit. However, using periodic boundary conditions instead 
of Dirichlet boundary conditions allows us to rely on a single periodic 
cell (Braides and Piat, 1995; Chambolle and Thouroude, 2009) via the 
following procedure. Let Y C R be the cuboid cell describing a single 
period of the periodic functions f and g. Using periodic boundary 
conditions, the effective crack energy may be computed via 


1 
hom : 
= f — ; 2.27. 
g (2; n) vESBV. periodic |Y | I glz, [v]; = vu) da ( ) 
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Notice that the two formulas for f*°™ and g*™ are fully decoupled, i.e., 
the formula (2.24) does not depend on g and (2.25) does not depend 
on f. This decoupling results from the different scaling of the surface 
and the bulk term within the functional (2.22), as well as the condition 
(2.21). If this condition is not fulfilled due to a certain scaling of f or g 
with the length scale n, a coupling of the two terms may be observed 
(Pellet et al., 2019). 

The work of Braides et al. (1996) on periodic homogenization has recently 
been extended to the case of stochastic homogenization by Cagnetti et al. 
(2019). This is of particular importance for the application to industrial 
materials, which often exhibit a randomness in their microstructure. 
Cagnetti et al. (2019) investigated the functional (2.22) with the same 
assumptions (2.21) but considered non-periodic, random and ergodic 
functions f and g. The parameter n in (2.22) serves as a length parameter 
of the microstructure indicating a scale separation in a random medium, 
see Fig. 2.9b for a schematic illustration. In their work, they showed that 
both the homogenized functional (2.23) and the cell formulas (2.24) and 
(2.25) are the same as in the periodic case. Hence, the decoupling of the 
two formulas holds upon stochastic homogenization as well. For the 
stochastic homogenization result the infinite volume limit for both the 
bulk term and the surface term is indispensable. 

Let us now consider the Francfort-Marigo model (Francfort and Marigo, 
1998) and discuss how the homogenization results of Braides et al. (1996) 
and Cagnetti et al. (2019) apply for this model. For heterogeneous C and 
y with characteristic length 7, the weak form of the Francfort-Marigo 
functional reads 


FM, (u) = ARE: (2) : Vuda | (2) dA. (2.28) 


For a given discretization in (pseudo-)time this functional is to be mini- 
mized in each time step with the constraint that the jump set $,, must 
contain the jump set of the previous time step. The weak solution of this 
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minimization problem is found in the function space SBD (Chambolle 
and Crismale, 2021). This functional is very similar to the functional 
(2.22). In particular, g(x, [v], n) = y(x) is a suitable option, since y is 
strictly positive. Furthermore, C is positive definite and thus a similar 
growth condition as (2.21) holds. However, we notice two deviations 
between the models. Firstly, the bulk energy of Braides et al. (1996) 
depends on the gradient of v, whereas in (2.28) the bulk term depends 
on the symmetrized gradient of u. Thus, instead of SBV in case of 
Braides et al. (1996), the function space SBD (Chambolle and Crismale, 
2021) is required. Only in the case of anti-plane shear these two terms 
coincide. Lifting this restriction, however, is subject of current research. 
In a recent preprint, Friedrich et al. (2022) investigated the case u € SBD 
for d = 2, i.e., lifting the restriction to anti-plane shear. For the 3D 
case they considered a homogeneous crack resistance. In their preprint 
they showed a periodic homogenization result of the Francfort-Marigo 
functional to the homogeneous functional 


FM" Vu: CË : Veude+ / yË (nu) dA. (2.29) 
Q Su 
The secend deviation between the Francfort-Marigo model and the func- 
tional (2.22) is given in the presence of a time stepping. The Francfort- 
Marigo model explicitly considers discrete time steps and introduces 
an irreversibility constraint. In the functional of Braides et al. (1996) 
on the other hand, no time stepping scheme is present. Without the 
irreversibility constraint, the results of Braides et al. (1996) and Friedrich 
et al. (2022) may be applied in each step separately providing a ho- 
mogenization result for the Francfort-Marigo model. This holds in 
particular for the initial time step if no pre-existing crack is present, 
i.e., Su is initially empty. The consideration of a time-stepping scheme 
and an irreversibility constraint in the homogenization result of Braides 
et al. (1996) has been provided by Giacomini and Ponsiglione (2006). 
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They showed a periodic homogenization result and found the same cell 
formulas as Braides et al. (1996). Therefore, in case of anti-plane shear, 
a periodic homogenization result for the Francfort-Marigo model is 
available. A proof for the general case, i.e., a stochastic homogenization 
result of the Francfort-Marigo model under anyy loading case including 
the irreversibility constraint, is, however, still pending, as the restrictions 
are only lifted individually. 

All aforementioned extensions of the homogenization result of Braides 
et al. (1996) show the same formulas for the effective properties and a 
decoupling upon homogenization. In order to compute the effective 
stiffness C°“ in case of periodic homogenization, one may rely on the cell 
formula (Bakhvalov and Panasenko, 1989), defined on a single periodic 
cell Y equipped with periodic boundary conditions. For a macroscopic 
strain field E € Sym(d) the cell formula reads 


E:C#:.E= inf es I (E + Vu): C(x): (E + Veu) dx. (2.30) 
u, periodic IY | Y 


Additionally, the case of coputing effective stiffnesses in stochastic 
homogenization is well-studied and any boundary conditions may be 
conducted. Upon an infinite-volume limit, the effects of the boundary 
conditions vanish (Sab, 1992; Bourgeat and Piatnitski, 2004; Owhadi, 
2003). However, for cells of finite size, the chosen type of boundary 
conditions does have an influence on the approximation quality of the 
"true" effective? stiffness, see the works (Sab, 1992; Kanit et al., 2003) 
among numerous others. It can be shown - both theoretically and 


5 A part of the mechanics community distinguishes apparent and effective properties. The 
former correspond to cells of finite size, whereas the latter emerge only upon an infinite 
volume limit (for stationary and ergodic media). Alternatively, apparent properties may 
be interpreted as approximations of the effective properties, in the same way as the 
displacement computed in a Galerkin discretization approximates the displacement 
of the continuous solution. In this work, we follow the second paradigm and use 
the terminology effective for quantities computed on cells of finite size, as well, tacitly 
assuming their approximative character. 
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numerically — that optimal convergence rates are reached when using 
periodic boundary conditions and periodized ensembles of microstruc- 
tures, see Schneider et al. (2022) for a thorough discussion. 

The cell formula for the effective crack energy is - in all aforemen- 
tioned cases — given by the formula (2.25). In the periodic case, the 
equivalent formula (2.27) applies. Thus, computing the effective crack 
energy breaks down to computing a y-weighted minimal surface, the 
minimum cut. 

Computing minimal surfaces in the context of fracture has been consid- 
ered for predicting crack propagation on two-dimensional micrographs 
even before the homogenization result of Braides et al. (1996) by Jeulin 
(1988; 1994a;b). In fact, in two spatial dimensions, the problem of com- 
puting the effective crack energy simplifies drastically. Indeed, it reduces 
to the problem of computing minimum (weighted) geodesics, for which 
efficient algorithms are available (Sethian, 1999; Osher and Fedkiw, 2002). 
Based on the homogenization result by Braides et al. (1996), Schneider 
(2020) established numerical tools to solve the cell formula for the 
effective crack energy on complex, three-dimensional microstructures. 


2.4.3 The cell formula for periodic minimum cut/maximum 
flow 


Consider the unit cell Y = [0, L1] x [0, L2] x [0, L3] and the field of crack 
resistances y : Y — R>o which satisfies for positive constants 71, ya > 0 
the condition 


yısyla)<Yy Vee y. 


Based on the homogenization result of Braides et al. (1996) and following 
contributions (Cagnetti et al., 2019; Friedrich et al., 2022; Giacomini 
and Ponsiglione, 2006) we define the effective crack energy via the cell 
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formula 


y(n) = inf, z) y(z)|n+Velde nes? (2.31) 
& periodic |Y] Y 
using periodic boundary conditions. This minimization problem seeks 
the periodic minimum cut through the cell Y with mean normal n. For a 
periodic cell Y this cell formula is equivalent to (2.25). A detailed discus- 
sion on the preferred boundary conditions in stochastic homogenization 
is given in chapter 6. 
A question of particular mathematical interest relates to whether a mini- 
mizer exists, i.e., if the infimum is actually attained, and if so, in which 
function space? Following Hintermüller et al. (2018) this is the case if 
7 is lower semicontinuous. We are mainly interested in heterogeneous 
materials for which y is a piecewise constant function and describes 
for instance inclusions of higher crack resistance within a matrix. To 
fulfill lower semicontinuity we set the value of y at the interface of the 
inclusions to the lower value. Under this assumption the minimizer is 
attained in BV. Notice, however, that the functional is convex but not 
strictly convex. Hence, any minimizer is a global minimizer but no 
uniqueness is guaranteed. Nevertheless, the resulting effective crack 
energy, i.e., the minimum value, is attained for any minimizer. 
The objective function of the minimization problem (2.31) is, for any 
normal n, given by 


1 
O7 I -y(ae) el) dz. 


This objective function is homogeneous of degree one in its argument, 
i.e, f(AS) = Af(€) VA > 0, and thus non-differentiable. Due to this non- 
differentiability of the objective function, gradient based methods to 
solve the minimization problem may not be applied directly. As a 
remedy we follow the strategy proposed by Schneider (2020). We 
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consider the constrained optimization problem 


1 f NE 
— | Ylig|| de — min with 
IY] Jy ECKE (2.32) 


Kp = {€:Y > B® |€=£+ V9, ¿e R?, ¢ periodic}. 


If the direction £ has length unity, the minimum value of this optimiza- 
tion problem is the effective crack energy for unit normal n = € given 
in (2.31). The problem (2.32) seeks the minimum cut € which satisfies 
the compatibility constraint specified via the set Kz. This constrained 
form permits a dualization of the problem, as suggested by (Schneider, 
2020). The formal dual problem to (2.32) is given by the maximum 
flow problem. This duality was first described by Strang (1983) who 
found that minimum cut is dual to maximum flow. The maximum flow 
problem seeks the periodic flow field v : Y — R?, solving 


1 _ 
— - v(x) de > max i 2.33 
IY] [é a div v=0, ||v(«)||<7(«) n 


This problem may be interpreted as a linear program with two con- 
straints. The first constraint is linear and enforces that the flow field is 
divergence free. This is the dual constraint to the compatibility of the 
cut field €. The second constrained bounds the local norm of the flow 
field v with the local crack resistance y. Both optimization problems, the 
minimum cut and the maximum flow problem are equivalent. Hence, for 
IE] = 1 both problems compute the effective crack energy with n = £. 

To solve (2.33) numerically, suitable discretizations and solvers are 
required. Schneider (2020) used an FFT-based solution framework 
with a trigonometric collocation discretization (Moulinec and Suquet, 
1994; 1998) and a finite element discretization with reduced integra- 
tion (Willot, 2015a). He solved the governing equations with a primal 
dual hybrid gradient method (Esser et al., 2010; Pock et al., 2009). A 
similar approach has been proposed by Willot (2020) who investigated 
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the related problem of finding the effective conductivity of resistor 
networks. Michel and Suquet (2022) proposed a different approach 
to computing the effective crack energy. For the discretization of the 
governing formulas they conducted a discretization method based on 
trigonometric collocation (Moulinec and Suquet, 1994; 1998). In their 
work (Michel and Suquet, 2022) they pointed out similarities of the 
governing equations with the problem of computing limit loads of 
structures (Christiansen, 1981). For their similar problem they conducted 
classical optimization methods. 

Finding novel discretization methods and solvers to compute the effec- 
tive crack energy via minimum cut/maximum flow is the main goal of 


chapter 4. 


2.4.4 Discussion on the terminology "effective crack re- 
sistance" 


Independent of the homogenization result for fracture, discussed in 
section 2.4.2 several approaches concerning effective properties in 
fracture mechanics are found in the literature. These methods define the 
effective crack resistance or the effective toughness using computations on 
heterogeneous materials and classical linear elastic fracture mechanics 
or phase-field fracture. Their approaches differ from our definition 
of the effective crack energy and we wish to put our approach into 
perspective. In particular, this difference gave reasons for naming the 
effective surface term when homogenizing the Francfort-Marigo model 
the effective crack energy instead of effective crack restistance. 

Bower and Ortiz (1991) provided a perturbative solution for a semi- 
infinite crack passing through a single, tough inclusion in a matrix. They 
relied on methods from classical linear elastic fracture mechanics and de- 
fined an effective toughness by either averaging or taking the maximum 
value of evaluated stress intensity factors during crack propagation. 
Roux et al. (2003) discussed an emerging effective crack resistance 
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for a material with isotropic and homogeneous elastic properties and 
a heterogeneous crack resistance. In this context, a self-consistent 
method for estimating the effective fracture toughness of a planar crack 
propagating through inclusions is established. For a medium with 
randomly distributed heterogeneities, they identified regions of weak 
pinning, where the fracture toughness is given by the arithmetic mean of 
the local toughness, and strong pinning, where a much higher toughness 
emerges, see also Démery et al. (2014) for a related study. Lebihain (2019) 
and Lebihain et al. (2021) extended the mentioned studies by accounting 
for cracks which bypass an inclusion, based on a perturbative, coplanar 
approach (Rice, 1985). 

To account for heterogeneity in the elastic properties, Hossain et al. 
(2014) performed phase-field fracture computations on heterogeneous 
microstructures with specific, so-called "surfing" boundary conditions. 
The emerging effective crack resistance equals the maximum in time of 
the J-integral evaluated along the crack tip, see also Kuhn and Miiller 
(2016) and Brach et al. (2019). 

Let us compare these approaches to compute the effective crack resis- 
tance with our approach based on the homogenization result (Braides 
et al., 1996; Giacomini and Ponsiglione, 2006; Cagnetti et al., 2019; 
Friedrich et al., 2022) for the Francfort-Marigo model of brittle frac- 
ture (Francfort and Marigo, 1998). The problem of heterogeneous 
fracture mechanics in general involves two prominent length scales: 
the correlation length of the heterogeneities and the typical size of a 
displacement increment. Typically, in a quasi-static framework, the 
size of the displacement increment is assumed to be infinitesimal. In 
this framework, when we consider a crack propagating through a 
microstructure, its progress may be hindered by various factors, like 
being pinned to an interface or avoiding an inclusion. This interpretation 
is implicit in Hossain et al. (2014) as well as Lebihain et al. (2021), for 
example, who consider crack propagation through microstructures of a 
fixed size and in continuous time, which is only discretized to enable a 
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numerical treatment. 

In practical applications however, the displacement increment is typi- 
cally of the order of magnitude of the macroscopic scale. In particular, 
in a homogeneous macro-model, the size of the heterogeneities has 
vanished due to a distinct scale separation and an incrementally 
discretized load is considered. This is the point of view considered 
when homogenizing the Francfort-Marigo model. The displacement 
increment is fixed, once and for all, as the limit of infinitesimally small 
heterogeneities is considered. Therefore, in this setting, the crack passes 
a microstructure cell within a single (macroscopic) increment. 

Another difference is the understanding of the emerging effective 
properties. As Hossain et al. (2014) always consider a time-continuous 
problem, their crack resistance is defined as the maximum in time of 
the local J-integral. Similarly, Lebihain et al. (2021) define the effective 
crack resistance by either the average in time or the maximum evaluated 
energy release rate. Notice that these definitions are neither related, nor 
motivated by mathematical homogenization results. Hence, in their case, 
a macroscopic material model has to be postulated. This macroscopic 
model then requires a discretization in time for numerical purposes. 

In contrast, the homogenization result for the Francfort-Marigo model 
practically works with an energy equivalence between the macro- 
scopic and the microscopic fracture energy, as a result of the energetic 
framework. In particular, no time step is present in evaluating the 
cell formula since the time step has been fixed once and for all on the 
macroscale. A detailed discussion on this difference is also found in 
Michel and Suquet (2022). 

Finally, let us remark that I-convergence implies the convergence 
of absolute minimizers, but does not predict what happens to local 
minimizers. Although an energy equivalence between the microscopic 
fracture energy and the macroscopic fracture energy appears in a 
natural way, the (absolutely) minimal surface in the cell problem is 
be a byproduct of I-convergence. From a physical point of view, it 
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might be more appropriate to work with crack surfaces that are just 
local minima of the weighted area. Still, keeping the assessment of the 
safety of microstructured components in mind, the absolutely minimal 
surface serves as a lower bound for the (real) effective crack energy, and 
is furthermore robust w.r.t. stochastic fluctuations in the microstructure. 
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Characterizing digital 
microstructures by the 
Minkowski-based quadratic 
normal tensor’ 


3.1 Introduction 


For material modeling of microstructured media, an accurate character- 
ization of the underlying microstructure is indispensable. The overall 
goal is to find simple quantities that describe the geometric shape as 
well as the composition of the microstructure under consideration. Once 
a microstructured material is characterized in these terms, correlations 
between microstructure characteristics and effective material parameters 
may be investigated via multi-scale methods, see Matouš et al. (2017). 

For microstructures of different material classes, different characteristics 
have been established. Particle reinforced composites may be characer- 
ized by the size(-distribution) of the particles and the volume fraction of 
the latter. The special case of fiber-reinforced composites additionally 
takes the fiber-orientation distribution, expressed via fiber orientation 


1 This chapter is based on Ernesti et al. (2022). In order to include this paper into the 
structure of this work I shortened the introduction and made minor changes to the 
manuscript. 
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tensors (Kanatani, 1984; Advani and Tucker, 1987) into account. For 
porous media, in addition to the porosity, i.e., the volume fraction of the 
porous space, the pore-size distribution (Kate and Gokhale, 2006), the 
turtosity (Neumann et al., 2019) and the chord-length distribution (Math- 
eron, 1975) are the typical quantities of interest. Polycrystalline materials 
are typically characterized in terms of a grain-size distribution (Dobrich 
et al., 2004) and an orientation distribution (Bunge, 1982; Böhlke, 2006). 
Since microstructures are often stochastic, a mathematical investigation 
of microstructure and their characterization is strongly related to the field 
of stochastic geometry. One basic tool in stochastic geometry to describe 
geometrical shapes and sizes is given by Minkowski functionals (Schnei- 
der and Weil, 2008; Hadwiger, 1951), also known as intrinsic volumes. 
They are defined for wide classes of shapes, including all convex sets 
and their finite unions as well as all bounded sets with smooth boundary. 
A Minkowski functional associates to any such shape a scalar quantity. If 
one requires such a functional to be invariant with respect to Euclidean 
motions, additive and to satisfy a certain continuity property, then it can 
be shown, see Hadwiger (1951), that in 3D it can be written as linear 
combination of only four basic functionals, the Minkowski functionals. 
Among them are the total volume, the total surface area, the Euler 
characteristic and one further functional, which for convex shapes may 
be interpreted as the mean width, or in the context of smooth boundaries 
as the integral of mean curvature. Approaches for computing Minkowski 
functionals are based, for instance, on marching squares (Mantz et al., 
2008) or on Steiner’s formula (Klenk et al., 2006; Guderlei et al., 2007). 

Being scalar-valued and rotation invariant, Minkowski-functionals are 
intrinsically insensitive to anisotropic features of the shape in question. 
Therefore, tensor-valued analogs of Minkowski functionals, the so- 
called Minkowski tensors (Alesker, 1999; Hug et al., 2008a;b; Jensen and 
Kiderlen, 2017), were introduced and studied. In addition to additivity 
and continuity, Minkowski tensors are required to be equivariant w.r.t. 
Euclidean transformations. This means, for instance, that rotating a 
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shape first and computing its Minkowski tensor afterwards leads to 
the same result as computing the Minkowski tensor first and rotating 
the tensor afterwards. A direct consequence of this property is that 
Minkowski tensors preserve axes of symmetry of structures, i.e., if a 
shape is rotationally invariant w.r.t. an axis p, the Minkowski tensor will 
be rotationally invariant w.r.t. p as well. 

Minkowski tensors may be computed for general microstructures with 
distinct interfaces, such as porous media, foams, bones or granular 
structures (Schréder-Turk et al., 2011). For porous media, Klatt et al. 
(2017) conducted a comparison between the common chord-length 
analysis and a Minkowski-tensor based approach. Schréder-Turk et al. 
(2011; 2013) evaluated Minkowski tensors for a given triangulation of the 
interface via explicitly known expressions for polytopes. Their ansatz 
was successfully used for characterizing the anisotropy of granular 
matter and metal foams, as well as identifying defects in molecular 
dynamics simulations of metal phases. For 3D gray-value images, Svane 
(2014; 2015) introduced approximation formulas for Minkowski function- 
als and tensors, also establishing convergence upon mesh refinement, 
called multigrid convergence in this context. Unfortunately, the cited 
works (Svane, 2014; 2015) did not include numerical examples. 

For finite point samples, Voronoi-based estimators Hug et al. (2017) may 
be used for approximating Minkowski tensors. 


Contributions 


We present an applied approach for characterizing digital microstruc- 
tures of industrial complexity in terms of the quadratic normal tensor, 
a tensor-valued quantity based on Minkowski tensors, bringing these 
concepts to the attention of the engineering community. 

For phenomenological continuum theories, which use microstructure 
information as state or microstructure variables to model the influence 
of microstructure on macroscopic material behavior, the Minkowski 
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tensors are promising quantities, because they are in principle observable 
and can be effectively calculated from three-dimensional image data. 
The Minkowski tensors complement, e.g., the already widely used 
fiber-orientation tensors (Advani and Tucker, 1987; Kanatani, 1984), 
which approximate the tangent distribution of the fiber centerline, and 
tensorial texture coefficients (Böhlke, 2006; Böhlke et al., 2010), which 
describe the distribution of crystal orientations. 

We introduce the relevant Minkowski functionals and tensors in section 
3.2 and isolate among them those suitable for microstructure char- 
acterization. In section 3.3, we present a novel algorithm for com- 
puting the quadratic normal tensor. For large microstructures with 
complex geometry, finding triangulations of the interface may be a 
challenging task, in particular if the microstructure is described by 
voxel data. Therefore, we present here an alternative to triangulation- 
based algorithms (Schréder-Turk et al., 2011) that works directly with 
gray-value images as input. The outward-pointing unit normals on the 
materials interface are approximated by finite-difference gradients of the 
discretized characteristic function. 

We investigate multigrid convergence of our approach by numerical 
studies in section 3.4. For fiber-reinforced composites, we compare the 
quadratic normal tensor to the more conventional fiber-orientation tensor 
of second order (Kanatani, 1984; Advani and Tucker, 1987). We compare 
the accuracy of our approach to the commonly used structure-tensor 
based algorithm (Krause et al., 2010) for computing fiber-orientation 
tensors. Last but not least, we study the anisotropy of sand grains and 
porous sand-binder aggregates based on the quadratic normal tensor. 
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3.2 Using Minkowski tensors for describing 
microstructures 


3.2.1 Minkowski tensors 


We briefly introduce Minkowski functionals and Minkowski tensors in 
a form suitable for our purposes and restrict to the 3D case. We refer 
to Schréder-Turk et al. (2013; 2011) or the lecture notes by Jensen and 
Kiderlen (2017) for the general case. 

Consider a solid body, by which we mean a bounded, not necessarily 
connected set K in R? with sufficiently regular boundary 0K. Here 
regularity can mean smoothness or convexity of some form. For our 
purposes it will be completely sufficient to assume that K is polyconvex, 
i.e., K can be represented as a finite union of (not necessarily disjoint) 
convex sets. To gain insight into the morphology of K, a shape index 
ọ associates to any such set K a scalar value. If one requires the shape 
index y to satisfy some natural basic properties, namely invariance 
with respect to rigid motions, additivity (meaning that (K U L) = 
(K) + (L) — (K N L) for solid bodies K, L) and a certain continuity 
(for convex sets, and w.r.t. the Hausdorff distance, see e.g. (Schneider 
and Weil, 2008, §12.3)), then it is a well-known fact due to Hadwiger 
(1951) that y may be represented as a linear combination of only four 
basic functionals Vo,..., V3, known as Minkowski functionals or intrinsic 
volumes. The Minkowski functionals encompass the volume V = V3, 
the surface area S = 2V3, and two further functionals, V; and Vo, which 
in special situations can be interpreted as the total mean curvature and 
the total Gaussian curvature of the body K. The latter is proportional to 
the Euler characteristic of K, i.e., the genus of the surface OK, which is a 
topological invariant. Volume and surface area are computed by 


V(K) = f dV and S(K)= | „3 (3.1) 
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Figure 3.1: Illustration of the e-parallel expansion K of the shape K C R3. (Ernesti et al., 
2022) 


If the boundary OK is sufficiently smooth, then the local mean curvature 
H and the Gaussian curvature G (i.e. the average and the product of the 
principal curvatures) are well-defined at each boundary point and the 
total curvatures may be computed via 

1 2 1 

V(K)=- HdS and Vu(K) = — G d5. (3.2) 

T Jax 4r Jax 
Such integral representations are also available for non-smooth bodies 
when one replaces OK by an integration over the normal bundle of 
K (Zähle, 1986). For practical computations, the additivity property is 
essential, allowing to decompose complex structures into simple convex 
pieces and to treat these pieces individually. For convex shapes, the 
Steiner formula provides another way to characterize the Minkowski 
functionals and another idea how to compute them. 
Consider, for a convex body K and e > 0, the -approximation 


K. = {x € R? : ||x—y]|| < £ for some y € K}, 
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see Fig. 3.1 for a schematic illustration. The Steiner formula (Schneider, 
2014) states that the volume of K+ is a polynomial in e, whose coefficients 
are (up to some normalization constants) the Minkowski functionals of 
K: 


V(K.) = V(K) +e8(K) + 1e2V;(K) + TEVR), (3.3) 


This allows to recover the Minkowski functionals of K by computing 
volumes of a number of e-approximations and inverting the above 
formula, see Klenk et al. (2006). One can also use the fact that the 
e-approximations K. are smooth even if K is not, allowing to determine 
the Minkowski functionals V; and Vo by means of the limit procedure 


Vi(K) = lim V (K) and Vo(K) = lim Vo(K.). 
e—0 e>0 


While these approximation results follow from the continuity of the 
Minkowski functionals, e-approximation properties of more general 
classes of sets are discussed in Rataj (2006). For more background to 
our informal discussion, we refer to Schröder-Turk et al. (2013) and the 


references therein. 


Since Minkowski functionals are, by definition, invariant w.r.t. Euclidean 
motions or change of frame, they are insensitive to directional and 
positional information. Hence, they are inappropriate for detecting 
anisotropies in a shape K. For this latter purpose and other applications, 
a more general theory of tensor-valued shape indices has been devel- 
oped, which are covariant w.r.t. Euclidean motions, see Schréder-Turk 
et al. (2013). In analogy to Hadwiger’s theorem (Hadwiger, 1951) and 
restricting to R? @sym R? = Roxy tensors, there are only six linearly 
independent shape indices (in addition to the Minkowski functionals 
multiplied by the identity), see Alesker (1999) and in particular (Hug 
et al., 2008b, §4). For a (convex) body K with sufficiently smooth 
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boundary, these may be expressed as 


ii 
WEK) = | xexav, ww =s | x@xdS, 
K aK 
1 1 
W3"(K) == H(x)x @xd8S, W3°(K) =; G(x)x 8 x dS, 
3 OK 3 öK 
1 
ww =5 | n& nds, W? (K)=- H(x)n @ nds. 
3 OK 3 öK 


(3.4) 


Here, x denotes the position vector of a point in K (or OK) and n 
stands for the field of outward-pointing unit-normal vectors on OK. For 
Minkowski tensors, Steiner-type formulas based on support measures 
have been established, see Schneider (2000). Note that some Minkowski 
functionals can be recovered from Minkowski tensors. For instance, the 
surface area is given by the formula S(K) = 3tr(W)?(K)). 
Based on these Minkowski tensors, Schröder-Turk et al. (2011) introduce 
the eigenvalue ratios 
minycEw) [Al 
BW) = ae (3.5) 
max eE(W) |A| 

as scalar measures of anisotropy. Here W stands for any of the six 
Minkowski tensors defined in (3.4) and E(W) is the set of eigenvalues 
of the symmetric matrix W. Clearly, 8(W) € [0,1]. For W = WZ”, 
W° and W?”, the matrix W(K) is positive semi-definite in general 
(and this is also true for the other Minkowski tensors if K is a convex 
body), implying that all eigenvalues of W are nonnegative. In this case, 
(W) = 1 if and only if all eigenvalues are equal, i.e., if the tensor is a 
multiple of the identity. Note that smaller values of 6(W) correspond to 
a higher degree of anisotropy. 
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3.2.2 Minkowski-tensor based microstructure characteri- 
zation 


Heterogeneous materials often exhibit random variations in their mi- 
crostructure, but a resulting deterministic material behavior (Torquato, 
2002). For characterizing microstructures, we are interested in singling 
out a small number of tensor-valued descriptors that may in turn be 
used as input for homogenization schemes, see Klusemann and Svend- 
sen (2010) for an overview. These microstructure identifiers should 
preferably exhibit certain natural properties: 


1. Respect for symmetries: We seek microstructure identifiers that pre- 
serve symmetry information. If a microstructure possesses some 
symmetry, then this is typically reflected in the macroscopic material 
behavior. Therefore, identifiers should capture such symmetry. 

2. Robustness: To be of practical use, small changes in the microstructure 
should only result in small changes in the descriptor. 

3. Translation invariance: For homogenization, statistical homogeneity 
is essential (Torquato, 2002). Thus, our identifiers should be invariant 
with respect to translations of the shape K. Furthermore, we want to 
explicitly include periodic structures, as periodic homogenization is 
often used for studying random microstructures (Kanit et al., 2003). 

4. Universal applicability: Minimal assumptions on the geometry of 
the structure allow for general application on a variety of different 


microstructures. 
In the light of these criteria, Minkowski tensors are promising candidates 
for microstructure characteristics. 


1. Their covariant tensorial nature reflects the anisotropy and direction 
dependence of the structure in question. 


2. They are robust due to their continuity properties w.r.t. the Hausdorff 
distance. For example, if a sequence of convex bodies K,, converges 
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to a convex body K, as n — oo, then all their Minkowski functionals 
and Minkowski tensors converge as well. Similar results hold e.g. 
if a polyconvex set K is approximated by its parallel sets K., see 
Schneider (2014) and Rataj (2006). There are also stability results 
showing Holder continuity with exponent at least 1/2, see Hug and 
Schneider (2015). 


3. If translation invariance is required, then beside the Minkowski 
functionals among the above mentioned Minkowski tensors precisely 
W? and W9? are suitable. As computing the curvature of interfaces 
of 3D voxel images is not straightforward (Lenoir, 1997; Monga et al., 
1991), we restrict in this article to the volume V, the surface area S 
and the Minkowski tensor W,’”. 

4. The Minkowski tensors are not restricted to specific shape as- 
sumptions on K. Indeed, only minimal assumptions on K are 
required (Schréder-Turk et al., 2011). For any practical application 
it is probably sufficient to note that any set (however complex) can 
be approximated arbitrarily well by a polyconvex set on which 
Minkowski tensors are defined. The tensor W,’” under consideration 
can in fact be defined under much weaker regularity assumptions, 
e.g. for sets with piecewise smooth boundaries. This flexibility 
distinguishes them from other approaches, where geometric priors 
are required for characterizing microstructures. For instance, for 
fiber-reinforced composites, fibers are often assumed to be (locally) 
cylindrical. Such geometrical priors run into problems for fibrous 
microstructure where the fibers deviate from their original cylindrical 
shape. For instance, during injection molding, fibers may be bent or 
twisted (Heinecke and Willberg, 2019). As they are independent of 
priors, Minkowski tensors may be suitable for characterizing fibers 
with distinct curvature. 

In the field of microstructure characterization, K is often the set union 

of a multitude of bodies, for instance inclusions within a surround- 
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n 


Figure 3.2: Body K with outward-pointing unit-normal field n. (Ernesti et al., 2022) 


ing matrix material. In this context, we are interested in a tensorial 
anisotropy-measure, which is stable w.r.t. an infinite-volume limit, where 
the number of inclusions tends to infinity. Thus, we normalize WP? to 
obtain the quadratic normal tensor (QNT) 


FR) 2 


which, for a single body or microstructure K with sufficiently smooth 
boundary, may be written in the form 


1 
ONT(K) = ay jones. 


where again n = n(x) is the field of normal vectors on OK, see Fig. 3.2. 
For a geometric interpretation of the QNT, observe that for any vector 
€ € R? the expression 


(n @ n)gé =n (n: &) 
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describes the orthogonal projection of € onto the line spanned by the 
normal direction n at x. In this sense, QNT(K) may be interpreted as 
an average over the normal projections computed w.r.t. the uniform 
probability measure concentrated on the surface IK. (Note that the 
resulting average matrix is still symmetric and positive definite but does 
not represent a projection anymore.) 

The ONT(K) admits an additional interpretation from a mechanical 
point of view. Suppose the structure K deforms with a homogeneous 
stress ø. Then, contracting the stress tensor with the QNT 


1 
str |, (as 


computes the mean normal stress on the surface OK. 


QNT(K) : o = 


By construction, the QNT is symmetric, positive semi-definite and has 
trace 1. In particular, QNT(K) admits an eigenvalue decomposition 
with real-valued, non-negative eigenvalues 1, Aa and A3, which sum to 
1. In case of a convex K, certain eigenvalue combinations can directly 
be interpreted in terms of the resulting shape of K: A; > Aa = A3, for 
instance, indicates a rather flat shape within the plane perpendicular to 
the eigenvector corresponding to Aı. For Aı = A2 > A3 we expect K to 
be a needle expanded in the direction of the eigenvector associated with 
A3, see also Appendix A.2, where the QNT is computed for a cylinder, 
and the sand grain experiments in section 3.4.4. 

Another advantage of Minkowski tensors is that they are locally defined 
and therefore locally computable. Complex polyconvex shapes can be 
cut into simple pieces and each piece can be treated separately. Then 
the additivity allows to recover the Minkowski tensor of the whole 
body from the Minkowski tensors of the pieces, allowing for efficient 
computation and parallelization. 
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3.3 Efficient implementation for 3D image 
data 


3.3.1 Algorithmic overview 


Consider a (periodic) heterogeneous two-phase material on the domain 
Y = [0,L.] x [0,L,] x [0,L;]. The microstructure of the material is 
described by its characteristic function x : Y — {0,1}, defining the 
two phases N and Q; via Qo = {x € Y : x(x) = 0} and Qı = Y\9o, 
respectively. Our aim is to describe phase Qı using the Minkowski 
functionals and tensors V(Q,), $(Q1), WI” (Q1) and QNT(,). 

Note that Qı is unknown in practice, only CT images of Qı can be 
observed. u-CT data is typically stored as 3D gray-value voxel data. We 
interpret the voxel data as a mapping xn : Yn — [0, 1] from the discrete 
set Yp, comprising the centers of a regular voxel grid with voxel length 
h, to the unit interval representing gray values. The gray value xr(y) 
associated to a point y € Y, stands for the volume fraction of (2; in the 
voxel centered at y. The relation between x and its discretization xn 
is demonstrated in Fig. 3.3. Fig. 3.3a shows the characteristic function 
x of a ball. Fig. 3.3b shows the non-discretized ball with the regular 
grid Yn in the background. In Fig. 3.3c, we see the discrete characteristic 
function xp of this ball as a gray-value image. Note that, in general, the 
input data x, does not allow to recover the phases Qo and Q; exactly 
as the interface is blurred. Only in the limit as h — 0 the correct 
characteristic function and, therefore, the correct sets are recovered. 
For determining w? 2 and S, in addition the normal directions are 
needed. In a weak sense, the unit normal n of the set Q; at a boundary 
point is recovered by n = —Vy, whereas -Vx = 0 away from the 
boundary. This statement may be formalized in terms of functions of 
bounded variation (Ambrosio et al., 2000). Therefore, we will compute 
the gradient numerically and establish formulas for Wp? and S based on 
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(a) Characteristic (b) x with background (c) gray-value image xn 
function x grid 


Figure 3.3: Characteristic function of a ball and its discrete representation by a gray-value 
image on a regular voxel grid. (Ernesti et al., 2022) 


volume averaging, see section 3.3.4. To improve the gradient estimation, 
a smoothing of the characteristic function x, is applied beforehand. 
The algorithm for computing the Minkowski quantities from a given 
voxel image is summarized in Alg. 1. First, we apply an image filter Fo 
to the characteristic function xp. Secondly, we estimate the outward- 
pointing normal vector by computing the gradient g from the resulting 
smoothed image 77. Finally, the desired quantities V, S, w? ? and QNT 
are estimated. 


Algorithm 1 Computation of Minkowski quantities 


Ti Fo * Xh > Blur image with image filter 
g(x) = VaTh (x) > Compute gradient 
Compute V by (3.7) 

Compute S by (3.8) 

Compute WP? by (3.9) 

Compute QNT by WP? /tr(W??) 

return (V, S, W®?, QNT) 
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3.3.2 Smoothing by image filters 


Due to the reconstruction procedure, u-CT scans often exhibit artifacts 
and impurities. Furthermore, binary voxel-based images do not allow 
reconstructing interfaces accurately (Mantz et al., 2008). To deal with 
these issues, we apply an image filter to the discrete characteristic func- 
tion (prior to computing the gradient). As different filters (and different 
choices of parameters) are available, we will also address choosing an 
appropriate filter. Applying the filter is realized by convolving the 
image with a specific filter kernel F. The filter parameter ø controls 
the width of filtering, and the result is the filtered image Z7, given as 
the convolution 


T; =F, * Xh- 


In our implementation, the convolution with F, is implemented via fast 
Fourier transform (FFT) (Cooley and Turkey, 1965), see Alg.2. Notice that 
in some cases the Fourier-transformed filter kernel may be computed 
efficiently without using the FFT. 


Algorithm 2 FFT-based filter application 


1: Xn + FFT (xn) > Transformation of the characteristic function 
2 Fyt FFT(F,) > Transformation of the filter kernel 
3: TS (E) — Xn lE NF (£) > Multiplication in Fourier space for all 


frequencies & 
4 Te + IFFT (Z? ) > Inverse transformation 


We shall consider a dimensionless filter parameter ø and scale it 
by the voxel length h. As filter kernels, we consider a Gaussian 


69 


3 Characterizing digital microstructures by the Minkowski-based QNT 


kernel (Bredies and Lorenz, 2018) 


CORNE eee ae |; 
Go(x) (ho)3(2m) 2 P ( — 


and the characteristic function of the unit ball, scaled to integrate 
to unity, 


B,(x) = | mr lei Sho, 
0 otherwise. 


In Fig. 3.4, the effect of filtering by a Gaussian and a ball kernel, re- 
spectively, is shown for a 1D laminate structure discretized with a 
voxel length of h = 2 um for three different filter parameters o. The 
red curve illustrates the impact of the Gaussian filter, whereas the blue 
line represents the ball-filtered image. In the Gaussian case, the resulting 
image is smooth across the laminate’s interface. However, for larger o, 
not only the interface is blurred, but no region of black or white remains. 
In fact, due to its global support, this even holds for small o. 

The impact of the ball filter is completely different. The piecewise 
constant indicator function with jumps at the interfaces is transformed 


into a piecewise linear function with slopes +4-. Therefore, when 
applying the ball filter to a structure with diameter larger than 2ho, 
some region with Z? = 1 will remain. 
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— Xh — Go * Xh — Bo * Xn 


0 0.2 0.8 1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 


0.4 0.6 . . 0.6 0.8 1 
zin mm zin mm zin mm 


(a) o = 0.055 (b) o = 0.15 (0 o = 0.2 == 
h h h 


Figure 3.4: Influence of filtering with different kernels and widths ø. (Ernesti et al., 2022) 


3.3.3 Approximating the surface normal by finite differ- 
ences 


Computing the Minkowski tensor WP? (Q4) requires determining the 
(unit) normal vector field n on the surface 0Q;. We approximate the 
normal field n by computing the gradient vector field 


8= Vilp 


of the filtered image Z} numerically. Notice that g is dependent on 
the voxel length h. At boundary points, we consider n = —g/||g]| 
as the outward pointing unit normal, provided g # 0. We briefly 
discuss the choice of the numerical gradient-approximation method. 
Finite-difference approximations are a simple way for approximating 
the gradient of a function given on a regular voxel grid numerically. 
Suppose a function f : Y > R, x > f(x) is given. We consider 
three finite-difference discretization schemes for the partial derivative in 
e;-direction (i = 1, 2, 3): 
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1. first-order approximation by forward differences, i.e., 


f(x + hei) — f(x) 


h Aj é 
o f(x) © h ’ 


2. first-order approximation by backward differences, i.e., 


f(x) = f(x = hei). 


al f(x) = ; Ä 


3. second-order approximation by central differences, i.e., 


F(x+he;) - f(x - hei) 


al f(x) © = 


Since our numerical experiments are performed on periodic structures, 
we treat the boundary in a periodic fashion. Under certain regularity 
assumptions on the function to differentiate, the first-order approxima- 
tions converge linearly in h to the exact gradient, whereas the second 
order approximation converges quadratically as h — 0, see Olver (2014). 
However, some further differences arise, which we demonstrate by 
example, see Fig. 3.5. 

Consider the filtered gray-value image of a ball, shown in Fig. 3.5a. When 
computing the gradient via central differences, the symmetry of the 
structure is recovered in the symmetry of the gradient field, since Z7 (x) 
and n(x) are evaluated at the same position, see Fig, 3.5d. However, if we 
compute the gradient via forward or backward differences, respectively, 
this will not be the case. The forward or backward partial derivatives in 


direction e; are not evaluated at x, but at x+ h/2e;, respectively, the faces 
of the cell. In particular, the partial derivatives in different directions 
will also be located on different faces. The resulting gradient fields are 
shown in Fig. 3.5c and Fig. 3.5b, respectively. Both appear deformed and 
uneven compared to the central-differences approach. Furthermore, a 
diagonal offset is noticeable. Numerical tests show that, in our present 
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(a) Filtered image Z7 (b) ||Vn Zz ||, computed via forward 
differences 
im 
- 
w 
= 
(c) |V aZ; ||, computed via backward (d) |V Z7 ||, computed via central 
differences differences 


Figure 3.5: Filtered gray-value image of a ball (a) and norm of the gradient computed via 
the three different finite-difference approximations (b)-(d). (Ernesti et al., 2022) 
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setting, the gradient approximation based on central differences is more 
accurate than the other approaches and will therefore be preferred, see 
also Section 3.4.2. 


3.3.4 Computing Minkowski tensors 


In this section, we propose formulas for computing the volume (fraction) 
of Qı, the surface area of 90, and the Minkowski tensor WP? (N1). Their 
accuracy and multigrid convergence will be investigated by numerical 
means in Section 3.4.2. 

The volume is approximated by quadrature, more precisely, by the 
trapezoidal rule, via 


VQ) = XO ah. (3.7) 


xEYn 


Motivated by results from geometric measure theory, see Giusti (1984) 
(Defintion 1.6, Theorem 1.24 and Definition 3.3) and Maggi (2012) (Propo- 
sition 12.20), we approximate the surface area via 


S(%) & SS lge), (3.8) 


xeY, 


where the gradient g(x) is computed by finite differences and || - || 
denotes the Euclidean norm. Due to the relation $(Q1) = 3tr(W??(Q1)), 
we approximate the Minkowski tensor w? 2 by 


We) E X sga) Te — 
i a lec) + e 


(3.9) 


where e > 0 is a small constant used to avoid division by zero. The ap- 
proximation of the quadratic normal tensor QNT is computed from the 
approximation of WP? by dividing by the trace, as in its definition (3.6). 
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3.4 Numerical examples 


3.4.1 Setup 


The algorithms 1 and 2 (as well as algorithm 3 discussed in Section 3.4.3 
below) were implemented in Python 3.7 with Cython (Behnel et al., 
2011) extensions. Critical operations were parallelized using OpenMP. 
For the eigenvalue decomposition of the structure tensor, discussed in 
Section 3.4.3 below, we rely on LAPACK (Anderson et al., 1999). The 
computations were performed on a desktop computer with a 6-core Intel 
i7 CPU and 32GB RAM. 


3.4.2 Parameter selection and multigrid convergence 


The proposed algorithm depends on several basic parameters including 
the grid size, the gray scale depth, the type and width of the applied 
filter and the choice of the gradient approximation, which we are free to 
choose in order to tune the algorithm. In this section, we investigate the 
influence of these parameters and propose suitable choices. 

In practical applications, the continuous range [0,1] of gray values is 
replaced by a set of discrete colors 


{0,1} if p = 1, 
en. ifp>2 


cP = 


of depth p > 1. In this context, we consider the discrete characteristic 
function as the mapping xn : Yn — CP. For p = 1, we voxelize the object 
under consideration in a binary manner, by colorizing a voxel if its center 
lies inside the object. For p > 1, we compute the binary image on the 
finer grid h’ = h/p and determine the gray-value of a voxel of size h as 
the mean value of its p? sub-voxels, resulting in a gray-value image of 
depth p. We investigate a ball Br of radius R > 0 for different p. The 
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Minkowski quantities of Br are known exactly and given by 


4 3 
V(Br) =, S(Bn) = 4a, 
An R? 1 
w2(Bp) = an Id, QNT(Br) = 3 ld, 


see Appendix A.1 for a derivation of the expressions for W®? and 
QNT. Hence these quantities can be compared to the corresponding 
numerically determined quantities V~, S~, W* = w? 27 and QNT~. 
For W = w? = and QNT, we define error measures by 


—- W” = NT(Br) - QNT*(B 
p= WER -W* (Be) na p= INTER) -ONT BR) | 
IW(Br)|| | QNT(Br)|| 
where ||- || denotes the Frobenius norm. Since W is connected to the 


surface area via S = 3tr(W), the error E is directly affected by an error 
in computing S. In contrast, this latter error does not necessarily affect 
E, as both QNT and ONT™ have trace 1. 

We investigate the influence of the different gradient approximations, 
filter kernels and filter widths o, as well as that of the depth p of the 
initial gray-value image, and we examine multigrid convergence as 
h — 0 numerically. The effect of the filters on the initial gray-value 
image, depending on the image depth, is exemplified in Fig. 3.6. In 
Fig. 3.6a, we see a slice through the characteristic function of a ball with 
a regular grid in the background. In Fig. 3.6b and 3.6c, we see slices of 
the discrete characteristic functions of the ball for depths 1 and 4. The 
center of the ball does not lie in the center of a voxel, but was chosen with 
a slight displacement, which results in a more uneven representation of 
the ball in the discrete images compared to Fig. 3.3. 

Apparently, images with a higher depth give rise to a more accurate 
representation of a ball than binary images do. Fig. 3.6d and 3.6e show 
the image after applying a ball-filter with o = 1.2. We see that the 
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difference between depth 1 and 4 has become smaller, but remains 
visible. The filtered binary image (p = 1) seems more uneven than the 
filtered gray-scale image (p = 4). 


(a) x with grid (b) xn for p = 1 (c) xn for p = 4 


(d) 2: forp=1 (e) 12 forp=4 


Figure 3.6: Characteristic function x, discrete characteristic function x» and filtered image 
Tj, ? for a single ball using depth 1 and 4. (Ernesti et al., 2022) 


First, we study the influence of the gray-value depth p of the initial voxel 
image for different spatial resolutions. The structure under consideration 
contains a single ball of diameter 16m in a box of edge length 241m, i.e., 
the material has a volume fraction of 15.5%. For this first study, we omit 
using a filter and rely on central differences for the gradient estimation. 
Fig. 3.7a shows the computed volume fraction vs. D/h, the diameter of 
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the ball per voxel length, for several gray-value depths p. The binary 
image, i.e., p = 1, exhibits the largest error and oscillates around the 
correct value. Only for a high resolution above D/h = 10, the error is 
within reasonable bounds. For a higher depth, the volume fraction is 
accurate even for the lowest resolution. 

The computed surface area vs. D/h is shown in Fig. 3.7b. Using the bi- 
nary image without any filter overestimates the surface area significantly 
and does not converge. For p > 2, we see that the error is reasonable for 
a resolution of 4 voxels per diameter and higher. For higher resolution 
and higher depth, the surface-area computation is rather accurate, but 
systematically overestimates the correct value by about 2% and does not 
converge. The error E of the Minkowski tensor is shown in Fig. 3.7c. 
For p > 2, it is below 6% even for the second-coarsest resolution of 
D/h = 4 and stays below 3% at higher resolutions. Finally, the quadratic 
normal tensor QNT is the one among the computed characteristics which 
is computed most accurately, see Fig. 3.7d. For p > 2, the error is 
below 5% for all spatial resolutions. Additionally, for all image depths, 
multigrid convergence is visible. This suggests that, to some degree, the 
error of computing the Minkowski tensor results from the mentioned 
overestimation of the surface area. Indeed, since ONT differs from w? 2 
by its trace and tr(W®?) = 9/3, the error of computing the surface area 
present in w? ? cancels out to some extent in QNT. 

In a second series of numerical experiments we repeated large parts of 
the above tests using first-order gradient approximations, as described 
in Section 3.3.3, instead of central differences. Compared to the latter, 
both first-order gradient approximations induce much larger errors, 
exceeding 20%. Therefore, we will restrict to central differences for the 
remainder of the article. 

Finally, we examine the influence of different filter kernels. Fig. 3.8 
shows the surface area as well as the two tensor-error measures vs. D/h 
for gray-value depth p = 1 (binary) on the left and p = 3 on the right. 
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— correct value — p=1 — p=2 —— p=3 — p=4 
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(c) Minkowski tensor error E (d) QNT error E 


Figure 3.7: Volume fraction, total surface area and tensor errors E and E for the unfiltered 
image, i.e.,o = 0. The gradient was computed via central differences. (Ernesti et al., 2022) 
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We consider the ball filter B, and the Gaussian filter G,, both with filter 
parameters o = 1.2 and o = 2, i.e., for a filter width slightly larger than 
a single voxel and a filter width of 2 voxels. 

In general, the errors for the gray-value image are smaller compared 
to the binary image. Focusing on the surface-area computation, i.e., 
Fig. 3.8a and Fig. 3.8b, we notice that applying no filter is actually most 
beneficial for a low spatial resolution. For p = 3, this even holds up 
to D/h = 10. For p = 1, the surface area is strongly overestimated for 
higher resolution. Even for p = 3, no multigrid convergence is achieved, 
if the filtering step is skipped. To achieve convergence, the ball filter 
with o = 1.2 is the most accurate. For p = 1, the ball filter with o = 1.2 
appears to be the best choice for resolutions up to D/h = 10. Above 
that threshold, the choice o = 2 exhibits the smallest error. Nevertheless, 
the ball filter with o = 1.2 serves as a good compromise. For both gray- 
image depths, applying the ball filter leads to better results than applying 
the Gaussian filter for computing the surface area of the structure. 
Investigating the error E, see Fig. 3.8c and Fig. 3.8d, permits us to draw 
similar conclusions. Fig. 3.8e and Fig. 3.8f show that the filter choice 
plays a subordinate role compared to the image depth for the error E. For 
the binary image, the Gaussian filter with o = 2 exhibits the lowest error, 
staying below the threshold of 2% for all resolutions. For gray-value 
images, however, the error of computing the quadratic normal tensor 
is below 3% for all resolutions and filters, which is accurate enough for 
most applications. 


3.4.3 A short-fiber reinforced composite 


Characterization of fiber-reinforced composites 

Short-fiber reinforced composites enjoy great popularity owing to their 
high (mass-)specific stiffness (Jones, 1998). The local fiber alignment is 
strongly dependent on the manufacturing process (Chung and Kwon, 
1995). The effective material behavior of short-fiber reinforced com- 
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Figure 3.8: Surface area and tensor errors E and E for depths p = 1 and p = 3 comparing 
the filter choice. (Ernesti et al., 2022) 
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posites is anisotropic, in general, and strongly dependent on the lo- 
cal fiber orientation. Each fiber is interpreted as a straight spherical 
cylinder of length L and diameter D, axis-aligned with unit vector 
p. Frequently used microstructure characteristics for fiber-reinforced 
composite materials are the volume fraction, the aspect ratio L/D and 
the fiber-orientation tensors of second or fourth order (Kanatani, 1984; 
Advani and Tucker, 1987). For fibers of equal length and equal diameter, 
the resulting fiber-orientation tensors (of order 2 and 4) of a structure 
with N fibers and orientation vectors pı,...,Pn are defined by 


N N 
A= WLP OP and A= HL POP SPP 

For varying fiber length and diameter, similar expressions have been 
proposed in Bay and Tucker III (1992) based on length- or volume- 
weighted averaging. 
For a gray-value „-CT image, the fiber-orientation tensors of second 
and fourth order may be computed by a variety of methods, see Pinter 
et al. (2018). A popular approach uses the structure tensor (Krause et al., 
2010), see Alg. 3. Alternatively, fibers may be segmented individually, 
see Hessman et al. (2019) for recent work. 
To gain insight into the relation between the fiber-orientation tensor 
A and the Minkowski tensor W,”, we compare their expressions for 
a single fiber of length L and diameter D, oriented in direction p, see 
Appendix A.2 for the detailed computation: 


A=p®p, (3.10) 
D? Ë 
we --pep+z(i-per)| and (3.11) 
1 È 
T= —— Z| Id- ' .12 
QN ma per+g(u por)| (3.12) 
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For microstructures containing N fibers, A is computed by averaging the 

single-fiber expression (3.10). w? 2 is computed by summing (3.11) over 

all fibers. The resulting quadratic normal tensor QNT may be computed 

as a surface-area weighted average of expression (3.12). 

The fiber-orientation tensor and the quadratic normal tensor need to be 

interpreted differently: 

e For a single fiber K, the fiber-orientation tensor of second order is 
a singular matrix (of rank 1) describing the projection onto the fiber 
axis. In contrast, the Minkowski tensor W? (K) of a single fiber K is 
a full rank matrix, which arises as a weighted sum of the orthogonal 
projection onto the fiber axis and the complementary projection onto 
the plane perpendicular to this axis. 

e For high aspect ratios, i.e., for L > D, for the QNT, the prefactor 
in front of the complementary projection is much larger than the 
other prefactor. 

e Using the fiber-orientation tensor as a descriptor of a microstructure 
rests upon specific assumptions that are often not satisfied for real 
structures. Typically, fibers are not of equal length, because they break 
during the manufacturing process (Inceoglu et al., 2011). Furthermore, 
the assumption that fibers are straight cylinders is not met in most of 
the cases, as longer fibers bend during manufacturing and therefore 
exhibit curvature (Heinecke and Willberg, 2019). In such situations, 
the structure-tensor based computation of the fiber-orientation tensor 
still gives some tensorial quantity as output. However, interpreting 
this result as a fiber-orientation tensor may not be justified. 

The Minkowski tensors, on the other hand, are not restricted to specific 
geometric assumptions such as particular shapes. Therefore, for 
structures containing curved fibers of different lengths or mixtures of 
fibers with other objects etc., W??? is still a geometrically well-defined 
quantity. As Minkowski tensors are integrals of locally computable 
quantities, see (3.4), they are even well-defined locally on any piece of 
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A QNT for 4 = 10 QNT for 5 = 25 QNT for 4 = 50 
100 0.048 0 0 0.05 0 0 0.012 0 0 
#1 000 0 0476 0 0 04875 0 0 0494 0 
000 0 00.476 0 0 0.4875 0 00.494 
0.79 0 0 0.1379 0 0 0.1219 0 0 0.131 0 0 
#2 0 019 0 0 0.3946 0 0 0.3996 0 0 04024 0 
0 0 0.02 0 0 0.4675 0 O 0.4785 0 00.484: 
0.49 0 0 0.266 0 0 0.2607 0 0 0.258 0 0 
#3 0 049 0 0 0266 0 0 0.2607 0 0 0.258 0 
0 0 0.02 0 0 0.468 0 0 0.4785 0 0 0.484 
06 0 0 0.219 0 0 021 0 0 0.205 0 0 
#4 0 03 0 0 0.348 0 00349 0 0 0349 0 
0.001 0 00.433 0 00.441 0 0 0.446 
0.33 0 0 0.333 0 0 0.333 0 0 0.333 0 0 
#5 0 033 0 0 0.335 0 0 0.335 0 0 0.335 0 
0 00.33 0 0 0.332 0 0 0.332 0 0 0.332 


Table 3.1: Comparison of fiber-orientation tensor A and quadratic normal tensor QNT for 
microstructures of different orientation and aspect ratio. 


a complex geometric structure. In contrast, the fiber orientation tensor 
is a non-local quantity intrinsically tied to cylindrical shapes. 


An overview of how the fiber-orientation tensor compares with the 
quadratic normal tensor for varying aspect ratios is given in Tab. 3.1. 
For this study, we generated 5 x 3 different microstructures, each con- 
taining 20% fibers of equal length and diameter, using the sequential 
addition and migration algorithm (Schneider, 2017). This algorithm 
draws fibers from an angular central Gaussian distributions on the 
two-dimensional sphere (Tyler, 1987). Indeed, the set of possible angular 
central Gaussian distributions may be parameterized by the second- 
order fiber-orientation tensors, see Montgomery-Smith et al. (2011). 
Across the microstructures we varied the orientation distribution (5 
different ones #1 — #5) and the aspect ratio (3 different choices: L/D = 
10,25 and 50). For convenience, all matrices are chosen to be diagonal 
w.r.t. the standard basis {eı, ea, e3}. 

Microstructure #1 is composed of aligned fibers in eı-direction. The 
second microstructure lies almost entirely within the eı — e2-plane, 
with preferred direction eı. The almost planar-isotropic case in the 
eı — e2-plane is realized via microstructure #3. A general anisotropic 
case with preferred direction eı and least preferred direction e3 is given 
in case #4. And, finally, microstructure #5 shows the isotropic case. For 
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the isotropic orientation (#5), all tensors are nearly equal. The QNT for 
the almost planar orientation (#3) results in one larger (corresponding 
to the normal vector of the plane) and two equal smaller eigenvalues, 
indicating no preferred direction within the plane. The QNT of structure 
#2 exhibits three different eigenvalues. The largest is equal to the largest 
eigenvalue of #3. Of the two smaller eigenvalues, the smallest indicates 
a preferred direction. The same interpretation holds for structure #4. For 
the uni-directional case (#1), the largest eigenvalue appears twice, which 
indicates a planar symmetry in both planes normal to the corresponding 
eigenvectors. By the smallest eigenvalue, again a preferred direction 
is indicated. 

In contrast to the fiber orientation tensor A, the quadratic normal tensor 
varies also with the aspect ratio L/D of the fibers. This may also be 
seen from the eigenvalue ratio 3 of QNT, see (3.5), listed in Tab. 3.2. 
This scalar measure of anisotropy is smallest in case of a unidirectional 
orientation distribution (#1) and almost 1 in the isotropic case #5. The 
degree of anisotropy is amplified for higher aspect ratios, which results 
in a lower £. 


B for & =10 | B for 4 = 25 | B for 4 =50 
#1 0.1003 0.0503 0.0250 
#2 0.2943 0.2544 0.2343 
#3 0.5690 0.5448 0.5327 
#4 0.5055 0.4751 0.4598 
#5 0.991 0.9903 0.9899 


Table 3.2: The degree of anisotropy of the different structures considered in Tab. 3.1 
measured by means of the eigenvalue ratio 6(QNT) of the quadratic normal tensor, see 
equation (3.5). Apparently, the anisotropy does not only depend on the fiber-orientation 
distribution, but is also sensitive to the aspect ratio L/D of the fibers 
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Sensitivity w.r.t. inter-fiber spacing 


A well-known challenge when computing fiber-orientation measures 
on w-CT scans is the sensitivity w.r.t. spatial resolution, as well as 
overlapping or touching fibers (Wirjadi et al., 2009). In the following 
study, we investigate the influence of the inter-fiber distance. Using 
the sequential addition and migration algorithm (Schneider, 2017), we 
generated structures with 20% fibers of aspect ratio 25, containing a 
total of 1336 inclusions. The fiber-orientation tensor was chosen almost 
planar isotropic with A = diag(0.49, 0.49, 0.02). The minimum distance 
between the fibers compared to their diameter can be chosen as an input 
for the microstructure generator. We generated 6 microstructures with 
minimum relative distance varying from 1% to 50% . Volumetric views 
and transverse slices of three of these structures are shown in Fig. 3.9. 
For 1% relative distance, several bundles of touching or almost touching 
fibers are visible, whereas, for 50%, each fiber is comfortably surrounded 
by matrix material. All structures were voxelized with gray-value 
depth p = 2 and for three spatial resolutions of D/h = 4,D/h = 8 
and D/h = 12, resulting in volume images with 256°, 512° and 768° 
voxels, respectively. 

For this data set, we compare the surface-area computation and the errors 
for the tensors W®? and QNT. As processing options, we compare no 
filter and the ball filter B, with filter parameter o = 1.2. The central- 
difference approximation is used for the gradient. Fig. 3.10a shows the 
computed total surface area vs. the minimum fiber distance relative to 
the diameter for the three spatial resolutions under consideration, using 
the ball filter with o = 1.2. We observe two trends. Firstly, the surface 
area is generally underestimated for all spatial resolutions. However, 
we clearly see multigrid convergence. Furthermore, the error is smaller 
for the larger minimum distance. This observation conforms to our 
expectations, as the surface area of touching or almost touching fibers 
is not computed accurately enough by a gradient-based approximation. 
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Figure 3.9: Fiber-reinforced composite containing 1336 fibers of equal length and varying 
inter-fiber spacing. The structures were generated synthetically using the sequential 
addition and migration algorithm (Schneider, 2017). (Ernesti et al., 2022) 
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Figure 3.10: total surface area and the errors E and E plotted vs. the minimum relative 
distance between fibers. (Ernesti et al., 2022) 
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In Fig. 3.10b, we see the results of the surface area computation without 
applying any filter. The errors are, in general, lower than in the case of 
o = 1.2. However, neither multigrid convergence, nor a convergence 
as the minimum fiber distance increases is observed. Fig. 3.10c and 
Fig. 3.10d show the error E for both filter choices. Again, the results 
reflect the relative error of the surface area estimation. Fig. 3.10e and 
Fig. 3.10f contain the errors of the quadratic normal tensor for both 
filter choices. For no filter application, the error is below 4%, and for 
o = 1.2, it is even below 2% for all spatial resolutions and minimum 
inter-fiber distances. No clear trend w.r.t. the inter-fiber spacing is visible. 
Hence, the quadratic normal tensor QNT may serve as a microstructure 
descriptor that is robust w.r.t. small inter-fiber spacing. 

We compare our approach with the well-established structure-tensor 
method, see Algorithm 3, which we implemented into our code. Com- 


Algorithm 3 Computing the fiber-orientation tensor via the structure- 
tensor method (Krause et al., 2010) 


1: LP + Fo *Xh > Blur image with image filter 

2: g(x) VaT; > Apply discrete gradient 

3: I(x) + g(x) 8 g(x) > Compute local tensor 

4.1, Furl > Blur local tensor with second filter 

5 {Ai (x), vi(x)} <— Eig(Z,,(x)) > Local eigenvalue decomposition 
(sorted, smallest first) 

6: A~ = Dovey, Vi(x) 8 vi (x) > Extract local orientation tensor 


7: return A~/tr(A*) 


puting the fiber-orientation tensor numerically via the structure-tensor 
algorithm requires applying a second filter F,, with filter parameter u to 
the tensor field n(x) © n(x) (component-wise, this tensor field denotes 
said ‘structure tensor’). Pinter et al. (2018) recommend that for the filter 
parameter for the second filter should be larger than for the first filter, 
which should be rather small. In our case, this is best recovered by 
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Figure 3.11: Error of the structure-tensor based fiber orientation tensor computation using 
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(d) Fiber-orientation tensor error EA, © = 
1.2, u =6 


different filter parameters for the first and second filter. (Ernesti et al., 2022) 
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choosing the ball filter with small filter parameter (i.e., o = 1.2) or no 
filter (i.e., o = 0) as the first filter. For the second filter, we choose a 
Gaussian kernel G, with u = 3 and u = 6. To evaluate the accuracy 
of the method, we introduce the fiber-orientation tensor error measure 
(similar to E) 


|A- A*| 


Bis 
|All 


where A” is the approximated fiber-orientation tensor computed by the 
structure tensor approach. 

Fig. 3.11 shows the error of this method for the four filter combinations 
o = 0,1.2; u = 3,6. The error is below 9% for all structures and 
resolutions. The first filter width o = 1.2 results in a lower error than for 
a = 0. This holds for all spatial resolutions and fiber-distance thresholds. 
For the second filter, however, the optimal choice depends on the spatial 
resolution. The two finer resolutions benefit from a larger second filter 
and even exhibit a larger error for the smaller p than for the coarse 
resolution. With respect to the relative minimum distance of fibers, no 
clear trend is visible. The error fluctuates between 1% and 8% for the 
different microstructures. The error of the quadratic normal tensor QNT, 
on the other hand, was below 2% for all spatial resolutions and hence 
provides a reliable option for characterizing fiber-reinforced composites. 
All computations were performed in a matter of minutes. 


3.4.4 Sand grains and sand-binder composites 


For manufacturing parts with complex geometry, casting is often the 
preferred choice (Rao, 2003). For casting, the mold enters a cavity of 
the specified shape. This cavity, in turn, is realized as a sand core, 
which has to be destroyed after the casting process. Such sand cores 
are composed of sand grains which are held together by an organic 
or inorganic binder. These constituents, their proportion and shape, 
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strongly influence the overall material behavior of the sand-binder 
aggregate (Schneider et al., 2018). Loosely speaking, if the strength of 
the aggregate is too low, the part will not survive the casting process. On 
the other hand, excessive strength may prevent the part to be extracted 
unscathed from the sand core. 

In this section, we compute the quadratic normal tensors of sand cores 
to study their anisotropy and to demonstrate the wide range of applica- 
bility of quadratic normal tensors. We consider six different sand-grain 
shapes which were obtained from fitting cleaned up and binarized 
w-CT scans (Schneider et al., 2018). The individual grains are shown in 
Fig. 3.12. 


ie u 
| 
| / \ 


(a) Grain #1 (b) Grain #2 


— — ae S 


(d) Grain #4 (e) Grain #5 (f) Grain #6 


Figure 3.12: Six different sand grains whose shapes are analyzed using QNT, see Tab. 3.3. 
(Ernesti et al., 2022) 
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Grain QNT B Grain ONT B 
0.2743 0.0518 0.0722 0.3669  —0.0428 0.0405 
#1 0.0518 0.2486 0.0179 0.4044 | #2 —0.0428 0.2949 —0.0204 0.6643 
0.0722 0.0179 0.4772 0.0405 —0.0204 0.3382 
0.382 0.0405 0.0055 0.2195 0.0407 0.0319 
#3 0.0405 0.2945 —0.0462 0.626 | #4 0.0407 0.3979 0.0936 0.4244 
0.0055 —0.0462 0.3235 0.0319 0.0936 0.3826 
0.28 0.0092 0.0013 0.3123 0.07 0.0428 
#5 0.0092 0.2875 0.07 0.5564 | #6 0.07 0.3202 —0.023 0.5795 
0.0013 0.07 0.4325 0.0428 —0.023 0.3674 
Table 3.3: Quadratic normal tensor QNT and eigenvalue ratio ß of the six grains in Fig. 3.12. 


These sand grains are non-convex and anisotropic. The computed 
quadratic normal tensors QNT are listed in Tab. 3.3. For the com- 
putation, we chose the ball filter B, with o = 1.2 voxels and used 
central differences for the gradient-approximation. In addition to the 
quadratic normal tensor, we quantify the degree of anisotropy by listing 
the eigenvalue ratios (3.5) of QNT. We observe that all sand grains have 
a distinct degree of anisotropy, varying between p = 0.4 and 8 = 0.63. 
To gain further insight into the anisotropy of the grains, we compute the 
eigenvalue decomposition of QNT for grain #1. The eigensystem reads 


—0.3233 —0.6666 
Aı = 0.5046, vı = | -0.1308 |; A2= 0.2914, v2 =] -0.6715 |; 

—0.9372 0.3236 

—0.6717 

de = 0.204, v3 = | 0.7294 

0.1299 


The largest eigenvalue indicates a somewhat disc-like shape within 
the plane normal to vı. The vectors va and v3 lie in that plane, the 
lower eigenvalue 3 indicates a slight extension in direction v3. For 
a better understanding, we point at Tab. 3.1, where fiber-reinforced 
composites are analyzed and a ‘translation’ to well-known orientation 
tensors is provided. 
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(a) Structure #1, containing 216 sand grains (b) Structure #2, containing 343 sand grains 


Figure 3.13: Sand core structures, containing 58.58% sand and 1.28% inorganic binder. 
The structures were generated by the mechanical contraction method (Schneider et al., 
2018). (Ernesti et al., 2022) 


These six sand grains of Fig. 3.12 were used for generating sand-binder 
composite microstructures, characteristic for casting applications, using 
the mechanical contraction method (Schneider et al., 2018). Two realiza- 
tions, containing 216 and 343 sand grains, are shown in Fig. 3.13. Both 
structures consist of 58.58% sand grains and 1.28% inorganic binder. In 
contrast to particle-filled composites, these microstructures involve an 
interpenetrating porous phase. 

On u-CT images of sand-binder composites, the binder phase cannot be 
distinguished from the sand phase, see Schneider et al. (2018). There- 
fore, we investigate how the presence of the binder phase affects the 
Minkowski tensors. We compare the quadratic normal tensor QNT of 
the sand-binder composite to the one with only sand grains for both 
structures in Fig. 3.13. We chose the ball filter B, with o = 1.2 voxels 
and central differences for the gradient approximation. For all structures, 
the resulting tensor QNT, the degree of anisotropy ( and the total 
surface area are listed in Tab. 3.4. The quadratic normal tensor is almost 
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Sand grains alone Sand-binder composite 
QNT 3 Simm] QNT B S{mm?] 
| 0.3294 —0.0002 —0.0025 ( 0.3292 —0.0015 —0.0006 ) 


Structure #1 —0.0002 0.3412 —0.0094 0.9329 51.74 —0.0015 0.3396 —0.0074 


—0.0025 —0.0094 0.3293 —0.0006 —0.0074 0.3312 


0.3218 —0.0004 —0.0011 0.3233 —0.001 —0.0014 
Structure #2 —0.0004 0.3366 —0.0028 0.9385 80.72 —0.001 0.3362  —0.0026 0.9451 73.96 


0.9485 47.26 


—0.0011 —0.0028 0.3415 —0.0014 —0.0026 0.3405 


Table 3.4: Quadratic normal tensor, degree of anisotropy and total surface area for grain 
structures #1 and #2 with and without binder. 


isotropic in all four cases. Removing the binder phase leads to slightly 
more anisotropic quadratic normal tensors compared to the sand-binder 
composite. However, the change is marginal. Without the binder, the 
surface area of every grain is fully exposed, which results in a 9.5% larger 
total surface area in case of structure #1 and 9.1% larger surface area in 
case of structure #2. In general, we see that, although the grains within 
both structures are highly anisotropic, the resulting microstructure as a 
whole is almost isotropic. Hence, mechanical contraction of anisotropic 
shapes results in an overall isotropic microstructure. This conforms to the 
results of Schneider et al. (2018), where elastic homogenization studies 
on similar structures were performed. An isotropic approximation of 
the effective stiffness tensor was shown to be accurate. 


3.5 Conclusion and Outlook 


In this study, we proposed using Minkowski tensors, a tensor-valued 
generalization of the scalar-valued Minkowski functionals, for the 
analysis of microstructures given implicitly on voxel images. Due to 
their tensorial nature, Minkowski tensors naturally contain information 
about the anisotropy of geometric structures and can be incorporated 
into continuum mechanical or other physical modeling approaches. 

We provide an efficient and compact algorithm for computing the 
Minkowski tensor Wf? and the resulting quadratic normal tensor 
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(QNT) from 3D gray-value image data. This algorithm is based on image 
filtering and a numerical gradient computation. We demonstrated 
the multigrid convergence of our algorithm on a single-ball structure. 
Central differences and a ball filter with low filter parameter turned out 
to be the most accurate for binary images. For gray-value images of low 
resolution, skipping the filtering step may be beneficial. Furthermore, 
we demonstrated that the quadratic normal tensor is rather insensitive to 
errors in the surface area computation, thus providing a robust measure 
of microstructure anisotropy. 

For fiber-reinforced composites, we compared characterizations based on 
the QNT to the well-established fiber-orientation tensors. We compared 
our approach to the common structure-tensor approach and demon- 
strated the accuracy and robustness of the quadratic normal tensor. 
Finally, we studied the QNT of sand-core microstructures. Its applica- 
bility to complex grain geometries demonstrates the versatility of the 
Minkowski-tensor approach. 

In future applications, further Minkowski tensors may be used for 
describing and characterizing a variety of microstructures, including 
curved fibers, fibers of different length and diameter, mixtures of several 
different shapes within a matrix, or polycrystalline structures. For a 
robust curvature-approximation technique based on voxel-image data, 
for instance, the curvature-dependent Minkowski tensor W,’* may be 
computed, providing additional information on the microstructure. 
The Minkowski tensors of the second rank may reflect only three types 
of material symmetries: isotropy, transverse isotropy and orthotropy. To 
detect finer material symmetries, working with higher-order Minkowski 
tensors is necessary. Mickel et al. (2013) suggested using irreducible 
Minkowski tensors for anisotropy characterization, a decomposition 
of the surface-normal density into those of some basic shapes in 
the spirit of Fourier analysis. This approach may also be beneficial 
for fiber-orientation analysis. Moreover, the concept of Minkowski 
maps (Klatt et al., 2012; Goring et al., 2013) may allow studying the local 
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differences of the fiber orientation across an inhomogeneous medium. 
Last but not least, Minkowski tensors may serve as input for further 
studies. Similar to fiber-orientation tensor based mean-field models (Ben- 
veniste, 1987; Kehrer et al., 2018), models based on Minkowski tensors 
may be developed. The quadratic normal tensor is able to provide 
insights for structures containing curved fibers and may serve as a tool 
for investigating their mechanical behavior. 
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Chapter 4 


An FFT-based method for 
computing the effective crack 
energy of a heterogeneous 
material on a combinatorially 
consistent grid’ 


4.1 Introduction 


Francfort and Marigo (1998) revisited Griffth’s original proposition 
(Griffith, 1921) in a quasi-static setting in order to include crack nu- 
cleation and crack branching, utilizing a variational formulation. More 
precisely, for a given body Q and after a discretization in pseudo-time, 
they seek the displacement u and the crack surface $ as minimizers of 
the Francfort-Marigo functional 


FM(u,8) = + 


1m Veu(x) : C(x) : Veu(a) dx +f y(z)dA (41) 


S 


1 This chapter is based on Ernesti and Schneider (2021). In order to include this paper into 
the structure of this work I shortened the introduction and made minor changes to the 
manuscript. 
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under the constraint of crack irreversibility, i.e., that the crack set S must 
contain the crack set of the previous time step. Here, V*u denotes 
the symmetrized gradient of the displacement field, i.e., the strain 
tensor field, C refers to the stiffness tensor and y denotes the crack 
resistance. Some care has to be taken with the formulation (4.1), as 
Griffith’s original proposal concerns only critical points of the functional 
(4.1) including local minima, local maxima and saddle points, whereas a 
rigorous mathematical treatment (Chambolle and Crismale, 2019) of the 
Francfort-Marigo model (4.1) appears to be limited to global minimizers. 
Please note that the formulation (4.1) accounts for heterogeneities in a 
natural way. 

A pertinent numerical approach to minimize the Francfort-Marigo func- 
tional (4.1) was introduced by Bourdin et al. (2000), and founded what is 
now referred to as phase-field fracture, see Wu et al. (2020) for a recent 
review. In close analogy to the Ambrosio-Tortorelli approximation (Am- 
brosio and Tortorelli, 1990) of the Mumford-Shah functional (Mumford 
and Shah, 1989), the phase field model approximates the crack surface 
S via a smeared interface of width l and introduces a damage variable 
d. Owing to their ability to nucleate cracks and to produce complex 
crack patterns, phase-field fracture models were subject to a flurry of 
activity (Ambati et al., 2015). In particular, strategies to account for 
material anisotropy in the phase-field framework (Prajapati et al., 2020; 
Teichtmeister et al., 2017; Schreiber et al., 2009) were proposed. As for 
linear elastic fracture mechanics of anisotropic media, such models may 
require elaborate and expensive experimental techniques to identify the 
material parameters. 

To alleviate this burden, multi-scale methods, in particular homogeniza- 
tion approaches, proved to be very effective for elastic and hardening- 
type inelastic material behavior. We refer to Matous et al. (2017) for a 
recent overview. Relevant for this chapter is a mathematical homog- 
enization result of Braides et al. (1996) for the Mumford-Shah func- 
tional (Mumford and Shah, 1989). In case of anti-plane shear and 
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neglecting the irreversibility constraint, this homogenization result is 
applicable to the Francfort-Marigo model of brittle fracture. More 
precisely, in a quasi-static setting and after a discretization in pseudo- 
time, we consider a fixed periodic microstructure (with non-degenerate 
stiffness and crack resistance). Following Braides et al. (1996) we identify 
the T-limit for vanishing period as the functional 


F Mer (u, S) = Veu(x) : Cog : Vru(e) da +f Yett(n) dA, (4.2) 
2 Jos s 
where n denotes the unit normal to the crack surface S. Here, the (possi- 
bly anisotropic) effective stiffness tensor Cept arises from the usual elastic 
homogenization formula based on the classical cell problem (Milton, 
2002). The integrand ye of the surface term is a function of a unit 
vector, and may be computed by a corrector problem involving the local 
crack resistances only. This corrector problem may be interpreted as 
finding the y-weighted minimal surface with average normal n cutting 
the microstructure (Schneider, 2020). 
In particular, the volumetric and the surface energies decouple upon 
homogenization, as a result of the different scalings of these terms in 
the model (4.1). Please note that this volume-surface decoupling is 
a consequence of the assumed non-degeneracy of the integrands. In 
case of degeneracy, an interaction of the two terms is not excluded, see 
Barchiesi et al. (2016) and Pellet et al. (2019). 
Recently, the homogenization statement was extended to the case of 
stationary and ergodic random materials (Cagnetti et al., 2019). Further- 
more, Friedrich et al. (2022) showed the homogenization result (4.2) for 
linear elasticity (without the restriction to anti-plane shear). Let us also 
highlight that the effective model (4.2) also emerges when homogenizing 
the Ambrosio-Tortorelli approximation of the Francfort-Marigo model, 
i.e., phase-field fracture models, see Bach et al. (2021). 
A method for computing the effective crack energy for three-dimensional 
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solids based on the homogenization result of Braides et al. (1996) and fol- 
lowing contributions was proposed by Schneider (2020). The approach 
is based on a convex reformulation of the minimum-cut problem (Strang, 
1983) in terms of maximum flow. More precisely, a primal-dual hybrid 
gradient method (Esser et al., 2010; Pock et al., 2009) was used, extending 
previous FFT-based computational homogenization methods for thermal 
conductivity and elasticity. 


Contributions 


This chapter is concerned with computing the effective crack energy 
via an appropriate cell formula (to be discussed in Section 4.2.1) corre- 
sponding to the mathematical homogenization results (Braides et al., 
1996; Cagnetti et al., 2019; Friedrich et al., 2022). Please note that the 
computed effective crack energy may differ from the effective crack resis- 
tance, depending on the underlying length scales and loading scenarios 
(see Section 2.1, item 5, in Schneider (2020) for a discussion). Still, the 
computed effective crack energy gives rise to a lower bound for the 
effective crack resistance, and may thus be used for assessing the safety 
of components made of such composites. 

Previous work (Schneider, 2020) provided a computational approach for 
computing the effective crack energy using an FFT-based primal-dual 
hybrid gradient solver and two discretization schemes, i.e., trigonometric 
collocation and the rotated staggered grid. The approach (Schneider, 
2020) has two shortcomings. First, the solution fields are characterized 
by ringing or checkerboard artifacts, depending on the discretization. 
Moreover, the solver does not permit reaching a high accuracy for 
complex three-dimensional microstructures. Although the first short- 
coming concerns the discretization and the second issue is related to 
the numerical resolution, both effects are actually related. Indeed, as 
computing the effective crack energy involves a pointwise constraint 
on a vector field, discretization-related artifacts may interfere with 
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solver performance. 

The cell problem for computing the effective crack energy is closely 
related to the minimum-cut problem put forward by Strang (1983) in 
his analysis of the continuous maximum flow problem, see Section 4.2.1. 
To be more precise, the maximum flow problem is rooted in graph 
theory and seeks a feasible flow through a flow network that obtains 
the maximum possible flow rate under capacity constraints (Ford and 
Fulkerson, 1956). With the help of duality theory for linear programming, 
it may be shown that the maximum flow equals the minimum capacity 
of a cut disconnecting the source and the sink (Ford and Fulkerson, 1956; 
Elias et al., 1956). Strang (1983) proposed a continuum generalization 
of the graph-theoretic maximum flow problem, involving an incom- 
pressible flow field subject to a (continuous) capacity constraint. Similar 
to the graph-theoretic version, he established a duality result which 
equates the maximum flow rate with the capacity of a minimum cut. 
However, the continuous version is no longer based on linear duality 
theory, and requires more sophisticated mathematical tools. Previous 
work (Schneider, 2020) realized that the cell problem corresponding to 
the mathematical homogenization results (Braides et al., 1996; Cagnetti 
et al., 2019; Friedrich et al., 2022) may be interpreted as a minimum-cut 
problem, where the microscopic (heterogeneous) crack resistance is 
regarded as a (spatially varying) capacity. Then, the maximum flow- 
minimum cut duality, valid on a discrete level, is invoked to construct a 
suitable primal-dual solver. 

For the graph-theoretic maximum flow problem, a variety of efficient 
solvers is available (Ford and Fulkerson, 1956; Dinic, 1970; Edmonds 
and Karp, 1972; Goldberg and Rao, 1998). Unfortunately, these solvers 
are unsuited for the continuous maximum flow problem. Indeed, work- 
ing with a graph-theoretic discretization leads to so-called metrication 
artifacts, which do not vanish upon mesh refinement (Kolmogorov and 
Zabin, 2004). More precisely, for the graph-theoretic version, the capacity 
constraints are associated to the edges, whereas, when representing 


103 


4 Computing the effective crack energy on a combinatorially consistent grid 


the continuous maximum flow problem in terms of a finite difference 
discretization, the capacity constraints are associated to the nodes of the 
graph. As a remedy, Couprie et al. (2011) introduced the combinatorial 
continuous maximum flow (CCMF) discretization, whose node-based 
capacity constraints account for all adjacent edges in a suitable way. The 
CCMF discretization, to be discussed in Section 4.2.2, naturally avoids 
metrication errors, and may be implemented into standard convex 
optimization solvers (Boyd and Vandenberghe, 2004). However, as 
the authors remark themselves: "In 3D, our CCMF implementation is 
suffering from memory limitations in the direct solver we used, limiting 
its performances." (Couprie et al., 2011, Sec. 4.4.3) 

Our contributions are threefold. For a start, we propose using the 
CCMF discretization for the effective crack energy associated to the 
mathematical homogenization results (Braides et al., 1996; Cagnetti et al., 
2019; Friedrich et al., 2022), see Section 4.2 for details. In this way, the 
artifacts of the previously used discretizations (Schneider, 2020) are fully 
eliminated. Our second contribution concerns a novel FFT-based solver 
for the maximum flow problem in the CCMF discretization on regular 
periodic grids, filling the gap mentioned by Couprie et al. (2011). The 
proposed solver, see Section 6.4, is based on a doubling of the degrees 
of freedom per cell, which makes the nonlocal capacity constraint for 
the CCMF discretization local, at the expense of additional constraints 
enforcing compatibility of the flow field across the faces of the voxel grid. 
As abyproduct, we arrive at an expression for the minimum-cut problem 
that is much simpler than in Couprie et al. (2011), see Section 4.3.1. Then, 
the alternating direction method of multipliers (ADMM), pioneered by 
Michel et al. (2000; 2001) in conjunction with FFT-based methods, is 
used to extract the minimum cut, see Section 4.3.2. Last but not least, 
we study recently proposed (Schneider, 2021b), adaptive strategies for 
choosing the penalty parameter in the ADMM. Finally, we demonstrate 
the capabilities of our approach in applications of industrial size, see 
Section 4.4. We find that an adaptive parameter-selection strategy is 
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Figure 4.1: Schematic of a potentially minimal crack traversing a two-phase microstructure 
for prescribed normal €. (Ernesti and Schneider, 2021) 


critical for high performance and high accuracy, improving upon the 
standard ADMM used by Willot (2020), who treats the closely related 
graph-based maximum flow problem (not the continuous one). 


4.2 The effective crack energy of a heteroge- 
neous material 


4.2.1 Cell formulas for the minimum cut and the maxi- 
mum flow 


Let us consider a cuboid cell Y = [0, Zi] x [0, L2] x [0, L3], on which 
a heterogeneous field of crack resistances y : Y — R is given. For 


mathematical reasons, we suppose that there are positive constants 7+, 
s.t. the inequalities 


y-<yla)<Yy+ holdforall wey. 
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We define the effective crack energy yer (Braides et al., 1996; Cagnetti 
et al., 2019; Schneider, 2020), a function on the unit sphere S? C R?, by 


aoa tl - Eu 
veal) = I y l|E+ Vol dx, Ees, (4.3) 


see Fig. 4.1, where |Y| = Lı L2 L3 denotes the volume of the cell and 
the infimum is evaluated over all smooth scalar fields 6 : Y + R 
which are periodic, together with all their derivatives. This formula 
computes the (periodic) minimum cut through the cell Y with mean 
normal ¿as a y-weighted minimal surface. The scalar ¢ plays the role 
of a characteristic function jumping across the cut. Please note that the 
integrand in the right hand side of equation (4.3) is a function which 
is homogeneous of degree one, i.e., it satisfies f(A€) = Af(£) for all 
vectors € € R and all scalars A > 0. This contrasts with thermal 
conductivity (Milton, 2002), where a homogeneity of degree two leads to 
a linear Euler-Lagrange equation associated to the variational problem. 
For the problem at hand (4.3), additional complications arise. For a 
start, due to the one-homogeneity, the functional in the definition (4.3) 
is not differentiable. In particular, the first-order necessary conditions 
are (strongly) non-linear. Furthermore, the one-homogeneity permits 
localization to appear for minimizers of the variational problem (4.3). This 
localization is not unwarranted, as such minimizers actually represent 
minimum cuts through the microstructure (Strang, 1983), weighted 
by the crack resistance, and enables computing the effective crack 
energy by the local crack resistance averaged over the minimum cut. In 
fact, the minimum cut needs not be unique. However, the computed 
effective crack energy is unique as a consequence of the convexity of the 
functional to be minimized. 

To circumvent the inherent lack of differentiability characterizing 
the functional (4.3), dual and primal-dual formulations may be ex- 
ploited (Chambolle and Pock, 2016). As an example, the (formal) 
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dual to the variational problem is given by the maximum flow prob- 
lem (Strang, 1983) 


[re max , (4.4) 
Y| Y div(v)=0 


loll <y 
where the maximum is evaluated over all smooth and solenoidal vector 
fields v : Y — R? which satisfy the pointwise constraint 


lula) < y(x) for (almost) all x € Y. (4.5) 


Due to the non-negativity of the terms involved, the latter condition may 
also be recast in the form 


ve)? < y(x)? for (almost) all x € Y. (4.6) 


The dual problem (4.4) maximizes the total flow in direction E through 
the microstructure under the point-wise constraints (4.5). The advantage 
of the dual formulation (4.4) over the primal formulation (4.3) is that it 
represents a smooth (in fact linear) optimization problem with linear 
and quadratic constraints, for which powerful solution methods are 
available (Boyd and Vandenberghe, 2004). However, some caution is 
advised, as the primal (4.3) and the dual problem (4.4) are strongly 
dual in the continuous setting only for a continuous crack resistance 
y (Strang, 1983). As soon as the crack resistance y is discontinuous, 
explicit counterexamples (Nozawa, 1994) to strong duality are known, 
i.e., the maximum computed in the dual problem (4.4) is strictly less 
than the minimum computed for the primal problem (4.3). 

For practical considerations, this delicacy does not play much of a role. 
Indeed, in finite dimensions, convex optimization problems with convex 
constraints always satisfy strong duality provided Slater’s condition is 
satisfied (Boyd and Vandenberghe, 2004, Sec. 5.2). Slater’s condition 
states that there is a strictly feasible point, i.e., a point where all inequality 
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Figure 4.2: Consistent placement of the flow-field variables on a generic voxel cell. (Ernesti 
and Schneider, 2021) 


constraints are satisfied as strict inequalities. Due to our prerequisite 
y > y- > 0, the field v = 0 is strictly feasible for the dual problem 
(4.4), and strong duality holds upon discretization. In particular, we 
may exploit the maximum flow formulation (4.4), as long as it arises by 
formal Lagrangian dualization (Boyd and Vandenberghe, 2004, Ch. 5) of 
a discretization of the cell problem (4.3). 


4.2.2 The combinatorial continuous maximum flow dis- 
cretization 


In this section, we discuss the combinatorial continuous maximum flow 
discretization (CCMF) introduced by Couprie et al. (2011) for the special 
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case of regular grids and in the periodic setting. The discretization 
scheme naturally approximates the continuous maximum flow formula- 
tion (4.4), and we take it as our point of departure. 

For this purpose, suppose that the unit cell Y = [0, L1] x [0, L2] x [0, Ls] 
is discretized by a regular grid with N; (i = 1,2,3) voxels for each 
coordinate direction. Each voxel is assumed to be cubic with edge length 
h, i.e., the conditions h = L;/N; (i = 1,2,3) are assumed to hold. Ina 
finite volume discretization, where each individual voxel serves as a 
control volume, the flow between adjacent cells is quantified by a flow 
field v, which is located at the voxel faces, see Fig. 4.2. The conservation 
of mass is encoded by the balance of in- and outflow 


0 = vit dik] — ofthat] + oft da] = vlii 


+ vfi,5,k+4] — v[i jk- 3], 
where we tacitly assume the integer indices i, j, k to satisfy 
0<i< Ni, 0<j<Na and O<k<N3, 


and the equation (4.7) should be interpreted in a periodic fashion. Let 
us denote by 


Wisk] = y ((i + 3h, (j + $)h, (k + 3A) 


the evaluations of the crack resistance y at the voxel centers, which we 
sample on a discrete grid Yy. Then, for the CCMF-discretization, the 
constraint (4.5) is approximated by the Nı Nə N; constraints 


2 yli,j,k]? > vļi+4,j,k]? + vfi-4,5;K] (4.8) 


+ vli.i-4,r]? + vfi,j,c+4 


Here, the constraint (4.8) is associated to each cell, and accounts for 
all six in- and outflow variables located on the corresponding adjacent 
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faces, see Fig. 4.2. Compared to the continuous formulation (4.6), which 
involves a vector of dimension three, twice the number of terms is 
considered. This is compensated by adding a factor two on the right 
hand side. 

Then, for prescribed average crack normal E e S°, the CCMF dis- 
cretization approximates the maximum flow problem (4.4) by the 


maximization problem 


1 =- = - 
N NN. > Eq vfit4,5,R] + Ey v[i,5+4,8] + Ez vfij k+] 
Pak (4.9) 


max . 
v satisfying (4.7) and (4.8) 


With FFT-based solution methods, to be discussed in Section 6.4, in 
mind, we transform the natural finite volume formulation into a more 
compact representation that is simpler to manipulate algebraically. For 
this purpose, we regard the flow field v as a vector field located at the 
voxel centers, with the identification 


Vali,d,k] = v[i+4,j,k], 
vylij,k] = vfi j+4,k], 


Yz[i,5,k] = vfij,k+4]. 
We also introduce a (backwards) divergence-type operator div via 


(div” v) [i,j,k] = Ur lé,d,4] — Veli—1,5,4] + Vyli,j,k] 


— Vyli,j-1,k] + Vz[i j,k] — Uz [i,5,k-1]. 


Then, the mass conservation (4.7) is satisfied precisely if div v = 0 
holds. To encode the constraint (4.8), we introduce the backwards shift 
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operator S, which operates as follows 


Uz [t—1,79,k] 
S(v) [i,k] = | vyli,j—-1,4] | - (4.10) 


Uzli,j,k—1] 
Then, the constraint (4.8) is equivalent to the condition 
vesel? + [Sotil < 2y, (4.11) 


expressed in terms of the Euclidean norm of the involved vectors. Last 
but not least, let us introduce the L? inner product on such vector fields 


1 


= HNN 2 (Velid kÜnlij k] + Vylid kylig k] + Valikli) 


i,j,k 
(4.12) 
with corresponding norm || - ||z2. With this notation at hand, we may 


express the maximization problem (4.13) in the compact form 


(€,v) 12 — max ; (4.13) 


„av v=0 
lull? +1] Se||? <27? 


where we regard € as a constant vector field and the norm constraint 
is enforced at every voxel. In the latter formulation, the similarities 
(and differences) to the continuous formulation (4.4) become apparent. 
Indeed, both the objective function and the divergence constraint are 
discretized in the natural way. The norm constraint, however, is replaced 
by a "non-local" constraint which involves neighboring values of the flow 
field, as well. Please note that this is a feature rather than a bug, as the 
flow-field variables are naturally located on the voxel faces, whereas the 
crack resistance is associated to the voxel center. Instead of interpolating 
the flow-field variables, the CCMF discretization averages the squares 
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of the flow fields. Such an approach has its merits, as will become clear 
in Section 4.4. 


4.3 An FFT-based solver for the CCMF dis- 
cretization 


4.3.1 The primal formulation for the CCMF discretiza- 
tion 


On a voxel grid, we consider the maximum flow problem (4.13) 


(€,v) 12 — max (4.14) 


div” v=0 
lol? +1] Sv||?<24? 


in the combinatorial continuous maximum flow (CCMF) discretization. 
With FFT-based resolution in mind, we compute the corresponding 
Lagrangian dual, i.e., the associated minimum cut problem. 

For later reference please notice that the adjoint of the backward shift op- 
erator S (5.14) w.r.t. the inner product (4.12) is given by the (periodized) 
forward shift operator 


Vg li+1,9,k] 
S*(v) 5,4] | Vylig+1.4] |- (4.15) 


Vzli,5,k+1] 


In particular, as backward and forward shifting are mutual inverses, the 
equation 5” $ = Id holds in terms of the identity operator Id. 

The shift operator is non-local, which makes the inequality constraint 
in the maximum flow problem (4.14) non-local, as well. With computa- 
tional resolution in mind, we seek a local formulation that relies upon 
a doubling of dimension. For this purpose, we introduce the linear 
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extension operator A, acting on vector fields v via 


1 v 
(Av) = Va | | ; (4.16) 


and producing a vector field with six scalar components per voxel. Then, 
the problem (4.14) may be expressed in the equivalent form 


(£v): — max, (4.17) 
div” v=0 
Avl <y 
where the factor two in front of the crack resistance (4.14) was transferred 
into the A-operator (4.16) and the norm in the constraint refers to the 
Euclidean norm of vectors with six components. For later reference, 
let us remark that the adjoint of the operator A (4.16) w.r.t. the six- 


component version of the L? inner product (4.12) is given by 


A* 
W2 v2 


= | Se (4.18) 


in terms of the backward shift operator (4.15). In particular, it holds 
1 
A* Av = z” + S* Sv) =v, 


i.e., A*A = Id and ||A|| = 1 in operator norm. Thus, the operator A 
is an isometric embedding, and the operator A* is a left inverse to the 
operator A. In turn, the operator AA* is the orthogonal projector onto 
the image of the operator A. To complete the necessary notation, we 
define the indicator function ur of a set T via 


TE 0 vweT 


+00, otherwise, 
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which permits encoding a constraint to the set T in terms of an 
objective function. 

With the necessary terminology at hand, we turn our attention to 
deriving the Lagrangian dual of the maximum flow problem (4.14) in 
the constrained form 


(£, v) 12 = Uf div v=0} (v) —Ec, (w) == wa (4.19) 


where we denote by C, the set 
Cy = {w: Yn > R®||wi.s,w]| < ylij forall i,j,k}. (4.20) 


The associated Lagrangian function reads 


L(v, w,€) = (& v) 12 — taiv- v=0} (V) = te, (w) — (€, Av + w)z2 (4.21) 


in terms of the Lagrangian multiplier field £ : Yy — R°. To evaluate the 
dual function 


() = sup L(v, w, 6), 


we rearrange the expression of the Lagrangian (4.21) 


(£) — sup(&, v) 12 ~~ Lfdiv- v=0} (v) Ze (£, Av) 12 +sup(€, w)L2 = CE (w) 
U —— w 
=(A*E,v) 72 
_ J) wa Lige Viik len, € € Kez, 
+00, otherwise, 


in terms of the set of compatible normal fields 


Kg = {£ : Yn > R°|thereissome $:Yy > R, 


E (4.22) 
st. A*E=E+Vto}. 
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This set may be interpreted as follows. Fields € : Yy — Rô associate to 
every voxel six scalar values. These values are assigned to each face of 
the voxel. Please note that any face of the voxel mesh is thus assigned 
with two values, one for each adjacent voxel. Up to a factor V2, A*E 
refers to the face-wise average of these two values. The set Kz contains 
all fields, whose averages A*€ are compatible in the sense that they arise 
as the sum of a constant vector (which is fixed beforehand) and the 
gradient of a scalar field (one scalar per voxel). Thus, up to face-wise 
averaging, the compatibility constraint is similar to thermal conductivity 
or elasticity. The final form 


1 
— i,j,k ijk i 4.23 
NNJN; > jkl llês ll — an (4.23) 


ofthe dual to the CCMF problem (4.14) is remarkably close to the original 
minimum cut formulation (4.3), cf. the more involved formulas in section 
2.3 in Couprie et al. (2011). 


4.3.2 An FFT-based ADMM solver 


To proceed, we rewrite the optimization problem (4.23) as an equivalent 
convex program that is amenable to operator-splitting approaches 


re) + ol) — min (4.24) 
in terms of the convex functions 


FOG) = ill) and IE = az DY Hot MEt 
i,j,k 
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The starting point of operator-splitting approaches is the rewriting of 
the unconstrained problem (5.21) in constrained form 


Fis) + gle) — = (4.25) 


For solving the problem (5.22), we utilize the alternating direction 
method of multipliers (ADMM) (Glowinski and Marrocco, 1975; Gabay 
and Mercier, 1976), which was pioneered in the context of FFT-based 
methods by Michel et al. (2000; 2001), and applied to non-smooth 
optimization by Willot (2020). For this purpose, we investigate the 
augmented Lagrangian function 


Ly(é,e,v) = FE) + gle) + (v, £- e)r2 +È lE- ellz; (426) 


involving a penalization factor p > 0 and the Lagrange multiplier v : 
Yn — Rê. The ADMM is based on the three-term recursion 


er — argmin, Lp(§, er, vr), 


etl — argmin, L,(E**!, e, v") 


yktl = yk ate poe _ er), 


(4.27) 


’ 


Let us investigate the first line more explicitly, 
en argmin, Lp(€, ef, v") 

2 

L2 


= argmin, f(E) + (de + 5 IE | 


2 
k 


| 
eae oe 
p 


= ie 2 ial ot oF 
752 


Thus, €**1 arises as the orthogonal projection of the point e" — v*/p onto 


the set Ke, 
: 1, 
ger = Pr; G _ av) . 
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Let us write down an explicit expression for the projection operator Pk; 
For given w : Yn > R®, we seek € : Yy — R, s.t. 


E=Px,(w), ie, = argmine cyl —w\|Z2 holds. (4.28) 


With the help of the orthogonal projector P = AA* and its orthogonal, 
complementary projector Q = Id —AA*, we may decompose the vector 
field € 

E=fp+ég with p= PE and 9 =Q6, (4.29) 


s.t., by orthogonality, 


Elli» = l£rlli» + allie (4.30) 


holds. Then, we may express the set X; in the form 


Kz ={€: Yn > RÔ | there are ¢ : Yn + Rand 7: Yn > RÔ, 


7 (4.31) 
s.t. € = A(E+ VB) + Qn}. 


To show equivalence of definitions, let us take an element € € Ke 
according to the former definition (4.22), which was characterized by 
the defining constraint 

A*E=E+VTO. 


Then, by definition of &p, we observe 
Ep = PE = AA*E = A(E+ V*9), 


and € is contained in the set (4.31). Conversely, applying A* to an 
element € of the set (4.31) yields 


Ke * E + * — + 
TISA TN nn d, 
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which shows that € lies in the original set (4.22). Returning to the 
projection problem (4.28) 


E= argmingex.liE — w||ĝ2, 


we decompose w = wp + wo in the same form (4.29). Due to the 
Pythagorean theorem (4.30), we observe 


é= argming ec || = wi: 
= argmin ex, lé — wells + Ilka - wol. 


= argmin,, ,||A(E + V*e) - wplli2 + Qn - wollte, 


where we inserted the definition (4.31) in the last line. We observe 
that the optimization problems for the variables ¢ and n decouple. The 
problem for Qn is particularly simple, and is solved by Qn = wo = Qu. 
Using that A is an isometric embedding, the problem for ¢ becomes 


é+ vro- Aw]. — min: 
The corresponding critical point satisfies 


div [E+ Vito A*w] =0 div Vt¢=div A*w. 


The latter equation may be solved formally to give 
= (div Vt)'div~ A*w, 


where } denotes the Moore-Penrose pseudo inverse. Reinserting the 
found expressions into the definition (4.31), we find 


E= AE + AVt(div” Vt)tdiv” A*w + (Id —AA*)w, 
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which we may also write in the more convenient form 


Px,(w) = AE + (Id-AA* + ATA*) w de 
with IT = V*(div V*)idiv-. 


The second line (5.24) can be rewritten using Moreau’s identity (Cham- 
bolle and Pock, 2016, Eq. (3.8)) in the form 


ge |p pe = Pe, yh pe") /p, 
where Pc, is the orthogonal projector 


(Po, (w)) War = Yeskl whe d/Mevted al wesi >, 
wfli,j,k], otherwise, 


onto the constraint set C, (4.20). Thus, we are led to the following scheme 
E+] = Ag — (Id —AA* + AT A*) (v* — per) ; 
p 


ehtl — [v® + pee Pe. (vt + a) Jp, (4.33) 


yet = yk £ p (ert _ ere, 


A recent study (Schneider, 2021b) highlighted the importance of utilizing 
a damping factor and choosing the penalty factor p adaptively. For this 
purpose, we consider the modified scheme 


ghti/2 = Ae x (Id -AA* + ATA*) (vë — o" ek), 
get] (1 — ö)ekt2 —{1—2é4)e*, (4.34) 


eft! = [v* + p" a = Pes (v* + pr eo) /p*, 


peti = yk a pe (en _ ertl), 
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with damping 6 € (0, 1) and adaptive penalty parameter p”. In general, 
the over-relaxation 6 = 1/4 is recommended (Monchiet and Bonnet, 
2012; Moulinec and Silva, 2014; Schneider, 2021b). Simple choices for the 
parameter p* are based on the Lorenz-Tran-Dinh scaling (Lorenz and 
Tran-Dinh, 2019) 
calle 
lle*llz2 


or the Barzilai-Borwein scaling (Xu et al., 2017) 


p (wkol eke 


= 
Tr, 


L2 


and an additional safeguard (Schneider, 2021b, Sec. 2.5). In our compu- 
tational experiments, the latter two schemes outperform, both, constant 
penalty parameter p and residual balancing (He et al., 2000). 

Last but not least, let us stress that the operator T has an explicit form in 
Fourier space, see (Willot et al., 2014, Eq. 18) 


4.4 Computational experiments 


4.4.1 Setup 


The algorithm (4.34) was integrated into an existing FFT-based compu- 
tational homogenization code for thermal conductivity (Dorn and 
Schneider, 2019), written in Python with Cython extensions (and 
OpenMP). In the context of small-strain inelasticity, the implementation 
of the ADMM (4.34) and the memory-efficient computation of the 
penalty factor is discussed in Schneider (2021b). In the same paper, the 


convergence criterion 


is. $ 1 
|fe u gk+3 


pe = fl Mey, (4.35) 
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(a) ||v||, CCMF (left), rotated staggered grid (center) and Moulinec-Suquet discretization 
(right) 


4. 
3. 
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1. 
0. 


(b) || ||, CCME (left), rotated staggered grid (center) and Moulinec-Suquet discretization 
(right) 


Figure 4.3: Flow v and normal € fields on a cross section through a 64° single-sphere 
microstructure for € = es, Yball = 10 Ymatrix and different discretizations. (Ernesti and 
Schneider, 2021) 


for prescribed tolerance tol, is identified as suitable. All computational 
experiments were run on a desktop computer with 32GB RAM and 
six 3.7GHz cores, and on a workstation with 512 GB RAM and two 
Intel Xeon(R) Gold 6146 processors (12 x 3.20 GHz), respectively. If not 
mentioned otherwise, we will use ADMM with damping factor ô = 0.25 
and the Barzilai-Borwein adaptive choice for the penalty factor. The 
default tolerance tol (5.33) was set to tol = 10%. 
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4.4.2 A single spherical inclusion 


As a first example, we build upon previous numerical experiments 
(Schneider, 2020, Sec. 4.2.2) and compare the CCMF-discretization to 
previously investigated discretization schemes, namely the rotated stag- 
gered grid (Saenger et al., 2000; Saenger and Bohlen, 2004; Willot, 2015a) 
and the Moulinec-Suquet discretization (Moulinec and Suquet, 1994; 
1998). We consider a 64° box containing a single spherical inclusion 
with a diameter of 32 voxels. The crack resistance of the inclusion 
is chosen as sphere = 10 Ymatrix- We prescribe a unit vector E = €y 
in x-direction as the crack normal. We solved the problem up to a 
tolerance of 10”? using ADMM and chose the penalty factor as lower 
bound p = min{7sphere; Ymatrix }, Which was the preferred choice for the 
primal-dual hybrid gradient method (Schneider, 2020, Sec. 3). Solution 
fields on a central cross section are shown in Fig. 4.3. 

The local flow fields v are shown in Fig. 4.3a. The Moulinec-Suquet dis- 
cretization shows significant artifacts, which is characteristic for Fourier 
spectral discretizations. The rotated staggered grid discretization, on the 
other hand, features checkerboard artifacts, although at a lower degree. 
In contrast, the flow field corresponding to the CCMF-discretization 
is much smoother, similar to the explicit jump discretization in the 
context of thermal conductivity (Wiegmann and Zemitis, 2006; Dorn and 
Schneider, 2019). The differences in the local crack-normal field £ for the 
CCME and the rotated staggered grid discretization are negligible and 
differ from the Moulinec-Suquet discretization in the local maximum 
values close to the central inclusion, see Fig. 4.3b. 

All discretization methods give rise to the same effective crack energy, 
i.e., Yet = matrix, aS the crack bypasses the inclusion in a plane. This is 
independent of the material contrast, as long as the crack resistance of 
the matrix exceeds the crack resistance of the single sphere (Schneider, 
2020, Sec. 4.2.2). 

As the Moulinec-Suquet discretization shows the strongest artifacts, we 
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focus on the remaining two discretization methods for the remaining 
investigations. 


4.4.3 A continuous-fibe reinforced composite 


In this section, we wish to assess the performance of the ADMM solver 
introduced in Section 6.4. As a measure of verification, we choose a 
comparatively simple microstructure which enables us to employ a high- 
fidelity interior-point solver (Domahidi et al., 2013) for second-order 
cone programs (Ernesti et al., 2021). The latter produces high-precision 
solutions, but is limited in terms of problem size. 

Accounting for this limitation, we consider a continuously fiber-rein- 
forced composite with 50% filler content. The two-dimensional mi- 
crostructure, containing 32 circular inclusions, was generated by the 
mechanical-contraction method (Williams and Philipse, 2003) and dis- 
cretized on a 128? voxel grid. The inclusions were furnished by a 
crack resistance of Yiber = 10 Ymatrix: We investigate the effective crack 
energy in direction € = e, and compare the CCMF discretization and 
the rotated staggered grid discretization, as well as different ADMM 
damping parameters ô, namely ô = 0.25 and 6 = 0.5. Furthermore, we 
investigate different selection strategies for the ADMM penalty-factor, 
the lower bound p = min{ ber, Ymatrix }, preferred in Schneider (2020) 
and the Barzilai-Borwein scaling (Xu et al., 2017), as well as the scaling 
by Lorenz and Tran-Dinh (2019) and residual balancing (He et al., 2000). 
As announced earlier, we compare the effective crack energy to solutions 
obtained by the high-fidelity solver ECOS (Domahidi et al., 2013) applied 
to conic reformulations of the minimum-cut problem (4.23) for the CCMF 
scheme and the discretization on a rotated staggered grid. We assess the 
solver quality in terms of the relative error 
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(a) ||v||, CCMF (left) and rotated staggered grid discretization (right) 
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(b) ||€||, CCMF (left) and rotated staggered grid discretization (right) 
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Figure 4.4: Cross section through the solution fields v and € for a 128? microstructure, 
containing 32 circular inclusions for CCMF and rotated staggered grid discretization for 
E = ex and fiber = 10Ymatrix. (Ernesti and Schneider, 2021) 
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in the effective crack energy, where Aurate is computed by the interior- 
8y: Teff P y 


point solver (Domahidi et al., 2013) with a residual of 107 +°. 

Fig. 4.4a shows the local flow field for, both the CCMF discretization and 
the rotated staggered grid discretization. For the rotated staggered grid, 
the flow field exhibits significant checkerboard artifacts in the inclusions 
as well as the matrix. The CCMF solution, on the other hand, is devoid 
of such artifacts. The corresponding crack paths are shown in Fig. 4.4b. 
The cracks bypass the inclusions and look qualitatively similar for both 
discretizations. However, the rotated staggered grid discretization 
shows a wider crack path, whereas the CCMF crack path is sharper. This 
allows the CCMF crack path to avoid several inclusions in a straight 
line, whereas the rotated staggered grid crack path has to avoid them, 
resulting in a less straight crack path. This observation is also reflected in 
the resulting effective crack energy, i.e. Yeff = 1.021 Ymatrix for the rotated 
staggered grid and Yes = 1.014 Ymatrix for the CCMF discretization. 

Fig. 4.5a shows the residual of the solver vs the iteration count for the two 
strategies for selecting the penalty factor, two damping factors and the 
two discretizations under consideration. During the first 1000 iterations, 
all solvers behave similarly, with a slight advantage for the choice 
6 = 0.25. After 2000 iterations, all solvers result in a residual below 
1078. For the CCMF discretization, the ADMM solver with ô = 0.25 and 
Barzilai-Borwein penalty-choice speeds up at 1600 iterations and reaches 
the required tolerance of 10~° shortly thereafter. For 6 = 0.5 a similar 
acceleration occurs after slightly more than 8000 iterations. Selecting 
the lower bound for the penalty factor p does not reach the required 
tolerance within 10000 iterations. 

For the rotated staggered grid discretization, the Barzilai-Borwein 
penalty-factor outperforms the constant choice, as well. For this 
discretization, the difference between the two damping factor choices is 
much smaller than for CCMF. 

The investigations are supplemented by Fig. 4.5b, which records the 
associated relative error (4.36). 
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(a) Residual vs iteration count, CCMF (left) and rotated staggered grid (right) 
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(b) Error vs iteration count, CCMF (left) and rotated staggered grid (right) 


Figure 4.5: Residual and error measure (4.36) for CCMF and rotated staggered grid 
discretizations, comparing different solver parameters. (Ernesti and Schneider, 2021) 
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Indeed, the relation between the residual (5.33) and the error in the 
quantity of interest (4.36) is not directly apparent. Indeed, we know 
that convergence of effective properties is implied by convergence of 
the fields. However, the quantitative relation between these may only 
be determined by comparison to a ground truth. For the CCMF dis- 
cretizations, the relative error (4.36) correlates with the residual rather 
well, reaching an accuracy below 10~* at convergence. In contrast, the 
solution for the rotated staggered grid leads to an error of only 0.5%, ie., 
hits a "stall". 

In addition to the mentioned penalty-factor choices, we studied two 
further (less competitive) approaches, namely residual balancing (He 
et al., 2000), which is often recommended in the literature, as well as an 
the approach suggested by Lorenz and Tran-Dinh (2019), which proved 
to be promising in small-strain micromechanics (Schneider, 2021b). To 
increase readability, the residual and the error (4.36) were moved to 
Fig. B.1b of Appendix B. 

With this validation at hand, we restrict to the CCMF discretization in 
combination with ADMM, damping factor 6 = 0.25 and the Barzilai- 
Borwein penalty-factor for the remainder of this chapter. 


4.4.4 A fiber-reinforced composite 


After the necessary verification steps, we turn our attention to problems 
with a higher degree of complexity. We consider a short-fiber reinforced 
composite with 18% filler content. The synthetic structure contains 376 
fibers with an aspect ratio (length/diameter) of 20, and was generated 
by the sequential addition and migration algorithm (Schneider, 2017). 
The prescribed fiber-orientation tensor of second order (Kanatani, 1984; 
Advani and Tucker, 1987) was diag(0.75, 0.19, 0.06), i.e., the fibers lie al- 
most exclusively in the x-y-plane with a strong preference in x-direction. 
The fibers are discretized with eight voxels per diameter, resulting in a 
volume image with 256° voxels, see Fig. 4.6a. Since the computations 
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Figure 4.6: Microstructure and effective crack energy for the fiber-reinforced composite. 
(Ernesti and Schneider, 2021) 


on such large structures are costly, we first investigate the influence of 
the tolerance entering the stopping criterion (5.33). For a configuration 
Yüber = 50 Ymatrix, We Computed the effective crack energy in direction 
€ = ez. After 1000, 2500, 5000, 7500, and 10000 iterations, we take a look 
at the corresponding residual and the computed effective crack energy, 
see Tab. 4.1. We observe that, after 1000 iterations we reach a residual of 
almost 107? with a relative deviation in effective crack energy about 2% 
compared to the prediction after 10000 iterations. After 2500 iterations, 
the relative error is below 1% with a residual at about 6 - 1074. For more 
than 5000 iterations, the effective crack energy does not change in the 
third significant digit, whereas the residual decreases only slowly. 

To complement these numbers, we take a look at a cross section through 
the computed crack surface at different iteration counts, see Fig. 4.7. 
We observe an influence of the solver accuracy on the solution field €. 
Indeed, after 1000 iterations, several distinct crack paths are present 
in the vicinity of the solution. These different cracks, however, come 
with different "intensities", as well. This ambiguity is reduced after 
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(a) 1000 iterations (b) 2500 iterations (c) 5000 iterations 


Figure 4.7: Cross section through crack surface for Yfiber = 50 Ymatrix at different ADMM 
iterations. (Ernesti and Schneider, 2021) 


# iteration | residual | Yer/Ymatrix 
1000 | 1.7 -1078 2.94 
2500 | 6.1- 1074 2.89 
5000 | 1.9- 1074 2.87 
7500 | 1.0 -1074 2.87 
10000 | 5.5: 1075 2.87 


Table 4.1: Residual and computed effective crack energy with normal £ = ez depending 
on the number of ADMM iterations for a fiber-reinforced composite, see Fig. 4.6a, with 
material parameters Yfiber = 50 Ymatrix 


2500 iterations. Only after 5000 iterations, the solver finds a unique 
crack surface. 

Please note that, in general, we do not expect the minimum-cut problem 
(4.3) to have a unique solution. Rather, for the problem at hand, a 
unique crack is formed and, at low levels of the residual, additional 
cracks appear, compare also Section 4.1.2 in Schneider (2020). These 
vanish, however, at high accuracy. To balance accuracy and ensuing 
computational costs, we fix the tolerance to 5 - 10-4. 

Next, we investigate the resulting crack surfaces and effective crack 
energies corresponding to different crack normals, see Fig. 4.8. In 
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‘fiber / “Ymatrix # iterations 
10 404 
20 | 2701 
30 | 2633 
40 | 2626 
50 | 2717 


Table 4.2: Number of ADMM iterations for varying material contrast Yaber/Ymatrix 


e,-direction, the effective crack energy is highest. This is caused by 
a preferred fiber orientation in this direction, forcing the crack surface 
to bypass the numerous inclusions. In e,-direction, see Fig. 4.8b, the 
crack surface looks roughly similar. However, the crack needs to avoid 
fewer fibers, resulting in a lower effective crack energy. In e.-direction, 
the crack surface is almost straight, see Fig. 4.8b, resulting in the lowest 
effective crack energy. 

Last but not least, we investigate the influence of the material contrast 
on the computed effective crack energy in = e,-direction, see Fig. 4.6b. 
This contrast is responsible for the allowed crack-inclusion interaction- 
mechanisms. Indeed, for high contrast, the inclusions can only be 
avoided, i.e., inclusion bypass is the only viable option. In general, 
the particles’ anisotropy (encoded by the aspect ratio and the fiber 
orientation for the example at hand) and the filler content determine 
the threshold in contrast where only inclusion bypass is permitted. For 
the example at hand, Fig. 4.6b reveals that this threshold is roughly 
at a material contrast of 40. In Tab. 4.2, the influence of the material 
contrast on the ADMM iteration count is listed. For a contrast of 10, the 
solver requires 404 iterations to reach the desired tolerance. Above a 
contrast of 20, the iteration count stabilizes at approximately 2700. For 
lower contrast, it may be energetically more favorable to cross some of 
the inclusions. For decreasing material contrast, this inclusion-crossing 
mechanism occurs more frequently. 
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(a) E = Er, Yet = 2.88 Ymatrix (b) € = ey, Yet = 1.67 Ymatrix (c) E = €z, Yeff = 1.20 Ymatrix 


Figure 4.8: Crack surfaces for the Cartesian normals, material contrast Yfiber /Ymatrix = 50 
and the fiber-reinforced composite, see Fig. 4.6a. (Ernesti and Schneider, 2021) 


4.4.5 Microstructures with a monodisperse pore distribu- 
tion 


As our next example, we consider microstructures with monodisperse, 
spherical pores and varying degrees of porosity. For a porosity between 5 
and 50%, we generated microstructures with 200 spheres by the mechan- 
ical contraction method (Williams and Philipse, 2003), see Fig. 4.9. All 
structures were discretized on a 256° voxel grid. The solid material has 
crack resistance y, and the spherical pores are furnished with a vanishing 
crack resistance, resulting in an infinite material contrast. Please note 
that in the previous study (Schneider, 2020, Sec. 4.1.1), the pores were 
furnished with a non-vanishing (yet small) crack resistance to ensure 
robust convergence of the utilized solution scheme. Such a restriction 
appears unnecessary for the improved solution method presented in 
this chapter. 

The effective crack energy in direction £ = e, and the required ADMM 
iterations are listed in Tab. 4.3. Following physical intuition, the effective 
crack energy decreases for increasing porosity. If the crack were straight, 
its crack energy would be proportional to the in-plane porosity of 
the crack plane. For a curved crack, there is a competition between 
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(a) 5% porosity (b) 25% porosity (c) 50% porosity 


Figure 4.9: Crack surface through microstructures with varying porosity. (Ernesti and 
Schneider, 2021) 


porosity in % | Yer/y | iterations 
5 0.871 4986 
25 0.535 3300 
40 0.412 3327 
50 0.305 2358 


Table 4.3: Influence of the porosity on the effective crack energy and solver performance 
for varying porosity 


"maximizing the porosity" and remaining as straight as possible, see 
Fig. 4.9. The iteration count appears to be uncorrelated with the porosity. 
As a remark, we found the iteration count to be strongly dependent on 
the specific realization of the microstructure, in general. Thus, we expect 
that no such correlation may be inferred from a single sample, but would 
require a more elaborate study. 


4.4.6 Sand-binder composite 


In our final example, we examine the microstructure of a sand-binder 
aggregate which is characteristic for inorganically bound sand cores 
used in casting applications. The synthetic structure was generated by a 
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(a) Microstructure 


| 


(d) Case #3 - porous inclusions 


(c) Case #2 - grain-matrix composite 


Figure 4.10: Bond sand-grain microstructure and Crack surfaces through the structure for 
three different combinations of crack resistances for matrix, inclusion and binder. (Ernesti 


and Schneider, 2021) 
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Ymatrix in MPa-um | Yerain in MPa-um | Ypinder in MPa-um 
#1 0 1 1 
#2 1 10 1 
#3 10 1 10 


Table 4.4: Material parameters for the three cases under consideration 


Ye in MPa-um | iterations 
#1 0.074 3204 
#2 1.133 1711 
#3 3.246 3971 


Table 4.5: Effective crack energies and iteration count for the three cases under considera- 
tion 


mechanical-contraction type method (Schneider et al., 2018; Ettemeyer 
et al., 2020), and is shown in Fig. 4.10a. The microstructure consists 
of three phases: The sand grains (58.6%), connected by a binder phase 
(1.3%), and a third phase (40.1%). In the physical applications, the latter 
phase represents the pore space. We wish to utilize the microstructure to 
get insights for a number of physical scenarios, and we will refer to the 
third phase more generally as the "matrix" for reasons that will become 
clear shortly. 

The crack resistances associated to the phases are denoted by Ygrain, Ybinder 
and ‘matrix, respectively. To investigate the effective crack energy and 
possible crack surfaces through the microstructure, we consider three 
different parameter scenarios, where the single phases model different 
physical scenarios. The governing parameters for the three cases under 
consideration are listed in Tab. 4.4. The resulting effective crack energies, 
as well as the required iteration counts are listed in Tab. 4.5. In parameter 
case #1, the crack resistance of the grains and the binder are equal, 
and the matrix material corresponds to a pore space. The resulting 
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crack surface is shown in Fig. 4.10b. We notice that the crack is fully 
contained in the binder phase. The effective crack energy is reduced to 
7.5% of the crack resistance which grain and binder share. The second 
parameter case models the structure as a matrix material with tougher 
sand-grain inclusions. The binder phase is treated as additional matrix 
material. Fig. 4.10c shows the crack surface avoiding the sand-grain 
shaped inclusions. The resulting effective crack energy of the composite 
is 1.133 matrix. The third case deals with the same contrast, i.e., the binder 
phase is once again treated as additional matrix. This time, however, the 
sand-grain shaped inclusions are weaker than the surrounding material. 
The effective crack energy is 32% of the matrix crack resistance. Fig. 4.10d 
shows the crack surface crossing several grains in order to avoid the 
matrix phase as much as possible. 


4.5 Conclusions 


In this chapter, we presented a powerful FFT-based solution method 
for computing the effective crack energy of industrial-scale composite 
microstructures. Based on a homogenization result for the Francfort- 
Marigo model of brittle fracture (Francfort and Marigo, 1998) in an 
anti-plane shear setting, see Braides et al. (1996), a cell formula for 
computing the effective crack energy was investigated. This cell formula 
may be interpreted as a minimum cut / maximum flow problem (Strang, 
1983), which finds various applications, for instance in graph networks 
and image segmentation. Following Couprie et al. (2011), we considered 
the CCMF discretization on regular voxel data and integrated it into an 
FFT-based computational homogenization framework. In comparison to 
traditional spectral and finite-difference discretizations, we found the 
CCME discretization to significantly reduce artifacts in the local fields. 

For solving the discretized equations, we investigated the alternating 
direction method of multipliers (ADMM) with various adaptive strate- 
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gies, and found a damping parameter 6 = 0.25 combined with the 
Barzilai-Borwein penalty-factor choice to be the most effective. We 
demonstrated the applicability of our approach to various large-scale 
problems, considering complex microstructures, as well as large or 
even infinite contrast in the local crack resistance. The presented 
framework was implemented into an existing homogenization code for 
thermal conductivity, and, although we ran some computations on a 
workstation, all presented computations could be done ona conventional 
desktop computer. 
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Chapter 5 


Computing the effective crack 
energy of heterogeneous and 
anisotropic microstructures via 
anisotropic minimal surfaces’ 


5.1 Introduction 


Francfort and Marigo (1998) reformulated Griffith’s theory of crack 
propagation (Griffith, 1921) as a variational problem, introducing the 
Francfort-Marigo model of fracture. For a fixed domain 2 and a fixed 
time discretization, they seek minimizers, i.e., the displacement field u 
and the crack surface S across which u may be discontinuous, of the 
energy functional 


1 
FM(u, S) == Vĉu : C(x) : Vu dz +f y(x) dA. (5.1) 
2 Jos S 
The energy consists of a volume-energy part accounting for elastic 
deformations, and a surface-energy part quantifying the crack-surface 
energy. Physical assumptions of their model include S‘ C S'*! at all time 


1 This chapter is based on Ernesti and Schneider (2022). In order to include this paper into 
the structure of this work I shortened the introduction and made minor changes to the 
manuscript. 
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steps t, i.e., the crack surface may only grow. Notice that their formula- 
tion is based on global minimization for reasons of mathematical well- 
posedness, whereas Griffith seeks local critical points. The Francfort- 
Marigo model naturally includes heterogeneous material properties, 
as both crack resistance y and the stiffness tensor C may depend on 
the position. Furthermore, distinct material anisotropy may already be 
expressed in terms of an anisotropic stiffness tensor. Additionally, as 
noted by the authors themselves (Francfort and Marigo, 1998, eq. (17)) 
the surface energy may account for anisotropy by replacing the isotropic 
crack resistance in the surface energy by an anisotropic term q(x, n(x)) 
depending on the unit normal n of the crack surface. 

The prevailing tool for treating the Francfort-Marigo model computa- 
tionally is the phase-field model of fracture. The method, introduced 
by Bourdin et al. (2000), is motivated by the Ambriosio-Tortorelli ap- 
proximation (Ambrosio and Tortorelli, 1990) of the Mumford-Shah func- 
tional (Mumford and Shah, 1989), used in image segmentation. The 
phase-field model involves a length-scale parameter l and seeks mini- 
mizers u as well as d, namely the displacement field and the damage 
variable, of the functional 

q(x) 7 


1 
PFi(u,d) = 5 fa — d)°V*u: C(x): Viu+ 2 T 


+ ivaj? dz. 
Q 2 


(5.2) 


Chambolle (2004) showed that, for L — 0, the phase-field model 
T—converges to the Francfort-Marigo model. Additionally, the phase- 
field model shows similarities to nonlocal damage models (Dimitrijevic 
and Hackl, 2008; Bažant, 1991) and may be treated as such as long as 
l > 0 is regarded as a material parameter (Kuhn, 2013). 

The phase-field fracture model is rather popular in the scientific 
community, including a variety of contributions over the last decade, 
see Wu et al. (2020) for a recent overview. Of particular interest to 
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this chapter are contributions that account for material anisotropy. 
Investigations on modeling anisotropic fracture using an anisotropic 
stiffness but isotropic crack energy were carried out (Gmati et al., 2020; 
Shanthraj et al., 2017). An approach to incorporate tension-compression 
anisotropy into the context of an anisotropic stiffness was suggested 
by van Dijk et al. (2020). Clayton and Knap (2014) introduced a 
geometrically nonlinear phase-field model with an anisotropic crack 
resistance, which takes, in case of small deformation elasticity, the form 
1 y(x) Ta? 


=(1- d)? Vu: C(x) : Vout —~|— + 1Vd- MPVd| de 


PF(u,d) =} alle 


i 
(5.3) 


with MP(z) = p(x) ® p(x) + B(x) (Id —p(x) ® p(x)) and a field of unit 
vectors p. They applied their model at small deformations to simulating 
cleavage in polycrystalline ceramics (Clayton and Knap, 2015). The posi- 
tive cleavage parameter ß is introduced to penalize crack propagation 
within the plane perpendicular to the unit vector p. Further extensions 
to incorporate crystal plasticity and ductile fracture were reported (Na 
and Sun, 2018; Bryant and Sun, 2018). The tensor MP permits to model 
either one weak plane, or one tough direction. Introducing a general 
symmetric and positive definite tensor M, up to three different crack 
resistances, i.e, the eigenvalues of M, in the three eigendirections may 
be prescribed. More general approaches were proposed using a multi 
phase-field setting, see Nguyen et al. (2017), or a higher order phase-field 
method, using fourth order tensors (Teichtmeister et al., 2017; Kakouris 
and Triantafyllou, 2019; Ma and Sun, 2020; Li et al., 2015), to study 
polycrystalline materials. Pillai et al. (2020) proposed an anisotropic 
cohesive phase-field model to simulate the behavior of anisotropic fiber 
structures. Incorporating weak interfaces via cohesive elements was 
proposed by Rezaei et al. (2021). 

To account for the influence of the microstructure of a material to its 
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macroscopic behavior, multi-scale methods may be used, see Matouš 
et al. (2017) for an overview. For hardening material behavior, those 
multi-scale approaches are well understood. Softening materials, in 
contrast, where distinct strain localization may occur, are more challeng- 
ing (Gitman et al., 2007). 

Braides et al. (1996) established a periodic homogenization result for the 
Mumford-Shah functional (Mumford and Shah, 1989) which is closely 
related to the Francfort-Marigo model of brittle fracture (Francfort and 
Marigo, 1998). Consider the Francfort-Marigo model with periodic 
material properties C and y with a fixed discretization in time, i.e., the 
quasi-static case. Furthermore, we neglect the irreversibility constraint 
and consider the loading case anti-plane shear. Subjected to these as- 
sumptions, the Francfort-Marigo model converges, as shown by Braides 
et al. (1996), to the effective functional 


FM*(u, 5) = Vou: CH: Voudz +f y(n) dA (5.4) 
s S 

with effective, possibly anisotropic stiffness tensor C° and effective 
crack energy q% (n). Furthermore, they showed that the two effective 
quantities decouple, i.e., the local stiffness tensor has no influence on the 
effective crack energy and the local crack resistance does not influence 
the effective stiffness. A key ingredient for the homogenization result of 
Braides and coworkers was a formulation on a single unknown, namely 
the displacement field which is permitted to be discontinuous across 
specific crack surfaces. In this way, the headache concerning the distinc- 
tion of a displacement field and the crack can be avoided. Indeed, the 
(possibly jumping) displacement field can be decomposed additively 
into a macroscopic and a microscopic part. 

The work of Braides et al. (1996) was further extended to stochastic 
homogenization by Cagnetti et al. (2019), ensuring representative vol- 
ume elements to exist for the Francfort-Marigo fracture model under 
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anti-plane shear. Recently, the restriction to anti-plane shear was lifted 
by Friedrich et al. (2022), so the homogenization result holds in general. 
To compute the effective quantities, Braides et al. (1996) provide specific 
formulas. Computing the effective stiffness reduces to classical linear 
elastic homogenization (Milton, 2002), whereas the effective crack energy 
eff 
ca 
y-weighted minimal surface cutting the microstructure. Note that the 


n) associated to a crack with normal n may be computed as a 


homogenization results are based on the notion of [—convergence. 
Thus, they only concern global minimizers of both the original and 
the effective functional. 

Schneider (2020) proposed a method for computing the effective crack 
energy for heterogeneous, locally isotropic microstructures using a 
reformulation of the cell formula of Braides et al. (1996) into a convex 
optimization problem. The transformation of the cell formula relies on 
Strang’s minimum cut/maximum flow duality (Strang, 1983). The es- 
tablished optimization problem was solved using FFT based algorithms. 
Recently, Ernesti and Schneider (2021) proposed a discretization on a 
combinatorially consistent grid, which improves the solution quality, 
introducing an associated adaptive ADMM solver, see also chapter 4. 


Contributions 


This chapter extends the approach presented in chapter 4 to account for 
a locally anisotropic crack resistance, enabling to treat matrix-inclusion 
composites with anisotropic phases or polycrystalline materials. In 
section 5.2.1, we introduce the cell formula for the effective anisotropic 
crack energy and anisotropic phases. We describe the anisotropy in terms 
of a tensor, similar to approaches applied in phase-field fracture (Clayton 
and Knap, 2014). We transform the anisotropic minimum cut problem 
into an anisotropic maximum flow problem in section 5.2.2, which gives 
rise to a convex optimization problem. 

Powerful solution methods for maximum flow problems arose for 
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maximum-flow problems on graphs (Ford and Fulkerson, 1956). Due to 
metrication artifacts, these are not directly applicable to the continuous 
maximum flow problem at hand (Kolmogorov and Zabin, 2004). This 
may be overcome by a combinatorial continuous maximum flow (CCMF) 
discretization (Couprie et al., 2011), recently applied to computing the 
effective crack energy of solids (Ernesti and Schneider, 2021), see also 
chapter 4. We propose a way to account for anisotropic crack resistances 
within the CCMF discretization in section 5.3.1. As the anisotropy 
of the crack resistance is described in terms of a (symmetric positive 
definite) second-order tensor, the maximum-flow formulation involves 
the projection onto an ellipsoid. We express the governing projection 
operator as the solution of a constrained optimization problem in section 
5.3.3. Finally, we apply our anisotropic minimum cut/maximum flow 
approach to brittle materials with a distinct anisotropy. In section 5.4.2 
we investigate polycrystalline brittle materials, such as ceramics with 
distinct cleavage in each grain. Last but not least, we investigate fracture 
of carbon fiber reinforced composites, where each fiber itself shows 
strong anisotropy of transversely isotropic kind. 


5.2 Minimum cut/maximum flow with 
anisotropic crack resistance 


5.2.1 Cell formula for an anisotropic minimum cut 


Consider a unit cell Y = [0, L1] x [0, L2] x [0, L5] and a given field of 
heterogeneous and direction-dependent crack resistances 


y:Y xR? oR (5.5) 
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We assume that, for any microscopic point x € Y, the association 
R? > > (a, €) 


defines a norm on the vector space R”, Le, it is one-homogeneous, 
satisfies the triangle inequality and is non-degenerate. Moreover, we 
assume that there are positive constants y_, 7+, s.t. the inequalities 


y-S7(@,€) <y} holdforall2eY and €€ R* with gl = 1. 
(5.6) 


Under these assumptions, we define the effective crack energy as a 
function on the unit sphere 9° via 


| = 
veal) = int |7 E+ Velo) de, 6.7) 


where the infimum is evaluated over all smooth periodic fields ¢. Upon 
one-homogeneous extension, the effective crack energy ert gives rise to 
a norm on R?, as well. 

In this article, we specialize the form of the microscopic crack resistances 
considered to those which describe an anisotropic Euclidean norm. More 
precisely, for a field G : Y — Sym(3) of symmetric, positive definite 
crack resistance tensors, we consider the local crack resistance field 


yz, 8) = || GK]. 


Then, equation (5.7) becomes 


gal = 
Yett(€) = inf z) || G(x) [E + Vo(z)] | dz. (5.8) 
$ |Y| Jy 
The isotropic case (Schneider, 2020; Ernesti and Schneider, 2021), also 
described in chapter 4, where we consider the isotropic crack resistance 


field y : Y > R, is recovered via G(x) = y(x) Id. 
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5.2.2 Anisotropic maximum flow formulation 


We define the energy functional 
©=— | Iesıa 69) 
= — x ; 
IY | Jy 


and, for fixed € € R3, the set of kinematic constraints 


Ke= {63 Y + R? periodic | £ =€+V¢ 


for some periodic ¢: Y > R}. 


With the energy functional and the kinematic constraints at hand, we call 


f(E) > min, (5.10) 


EEKE 


the anisotropic minimum cut problem. For fixed normal £, the minimum 
effective crack energy (5.8) computes as the minimum value of this 
minimization problem. 
Treating this problem numerically is challenging, since the energy func- 
tional is homogeneous of degree one and thus non-differentiable. Ex- 
tending the isotropic case (Schneider, 2020; Ernesti and Schneider, 2021), 
see also chapter 4, we introduce a dual formulation, i.e., a corresponding 
anisotropic maximum flow problem (Strang, 1983). The formal dual 
problem to the minimization problem (5.10) is given by 

Fw) 


-gE vdr min, (5.11) 


where f* denotes the Legendre-Fenchel dual of f, given by 


eee ei sine i= 
Flo) = sup jpj |,E vde- SE = sup yy f €v- Gelde 6.12) 
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and the minimum (5.11) is evaluated over all periodic solenoidal fields v. 
Since the tensor G is symmetric and positive definite, and thus invertible, 
we transform the Legendre-Fenchel dual (5.12) 


4 _ 1 
f (o) = sup p |, €- v — IG El de 


1 o en 
= up wf E Gv- [glia 
e-cel¥| Jy 


o Ge s1, 


+oo otherwise. 


Thus, the Legendre-Fenchel dual f* equals the indicator function 


0 ve, 
tc = 
+00 otherwise, 


of the closed set 


C={v:Y > R, periodic | || G(x)“ v(2)|| <1 (5.13) 
for almost all x € Y }. l 


Since G is symmetric and positive definite, the set C is a convex domain. 
More precisely, the constraint in the definition of the set C describes an 
ellipsoid centered at the origin, see Fig. 5.1. Combining (5.12) with this 
expression for C, with arrive at the anisotropic maximum flow problem 


1 = 
— -vdx — . 
IY] n ~~ divo=0, ja: v||<1 
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1Y 


Figure 5.1: Schematic sketch of the ellipsoidal domain Cz, i.e., all vectors satisfying 
||(2)-t (a) | < 1, with eigensystem {v;} and semi-axis 7;. (Ernesti and Schneider, 
2022) 


5.3 Numerical treatment 


Ernesti and Schneider (2021) proposed to solve the maximum flow 
problem in the combinatorial continuous maximum flow (CCMF) dis- 
cretization (Couprie et al., 2011) by FFT-based methods (Schneider, 
2021a). The formulation is based on doubling the degrees of freedom. 
We present an extension of this strategy to account for an anisotropic 
crack resistance expressed via a heterogeneous field of crack resistance 
tensors G : Y + Sym(3). We refer to Ernesti and Schneider (2021), as 
well as chapter 4, as a general reference throughout this section. 


5.3.1 The CCMF discretization for anisotropic maximum 
flow 


We discretize the unit cell Y = [0, Lı] x [0, La] x [0, L3] on a regular cubic 
voxel grid Yy with N;, (i = 1,2,3) voxels in each direction under the 
assumption that each voxel is cubic with edge length h = L;/N;, (i=1,2,3). 
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We evaluate the crack resistance tensor at the voxel center, i.e., 
G [ijk] = G ((i+ $)h, (j + 3)h, (k + 5)h) 
and the flow field v : Y — R? on the faces of each voxel 
Uz [i,j,k] = Uz (ih, (j+ s)h, (k+ 3 
2 


i+4)h, jh, (k + 


Uyli,j,k] = Vy ( 


( 
Uz[t,9,k] = Vz ((i 


This placement of the flow field enables to encode the conservation of 
mass via the discrete backwards divergence operator div” 


(div” v) [ik] = Veli,j,k]— Ve [i-1,9,k] + dyli,s,k]—Vy[l,5-1,kl+Vz[i,5,k]—Vz[i,5,k-1]. 


To enforce the constraint ||G~' v|| < 1 in the context of the CCMF 
discretization, we introduce the nonlocal shift operator S 


Vy [i-1,9,k] 
S(v)fi,5,4] | vylis-LA] |- (5.14) 


Uzli,j,k—1] 
The constraint is then enforced via the inequality 
2 2 = 2 
ICs. "vr + [|G iia So isl] < 2. (5.15) 


To integrate the CCMF discretization of the anisotropic maximum flow 
problem into an FFT-based homogenization solver for heat conductivity, 
we introduce the extension operator A, given by 


° | (5.16) 
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With this notation at hand, accompanied by the inner product 


1 


(v, w) = —_ 
Ni NoN3 


N deli. kwalijk] + Vylij.klwylij,k] + Vwlij,k], Veli9.Al 
i,j,k 


we may rewrite the discrete maximum flow problem as 


(€,v) — max with G= = : (5.17) 
div” v=0, 0 G 
|Č A(v) || <1 


For the CCMF discretization, the set C from equation (5.13) reads 
Ce= {w : Yy > RÊ, | | Grn wise <1 for alli, j, kb . (5.18) 


The derivation of the associated discrete primal problem follows the 
same steps as the isotropic case, described in section 4.3.1. The discrete 
minimum cut problem with a tensorial crack resistance is given by 


1 a 
MNN SORIRE 5.19 
Ni N2 N3 2 Ic: IwIE| ain — nin (5.19) 


with set of discretely compatible fields 


Kz = {£ : Yn > R®|there is some 6: Yy > R, 


(5.20) 
st. A*E = E+ VIO}. 


The operator A* denotes the left inverse of A, i.e., A* A = Id holds, and 
V* refers to the discrete forward gradient operator. 


5.3.2 The alternating direction method of multipliers 


To solve equation (5.19) with the alternating direction method of multipli- 
ers (ADMM), we rewrite the discrete minimum cut problem equivalently 
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as 
FE Paley > (5.21) 


with the two convex functions 


1 
F(§) = ır.(E) and g(f) = NiN2N3 3 [|G Wr Elisi], 
i,j,k 


where ix; is the indicator function of the set Kz described in (5.20). 
The operator-splitting approach starts by rewriting the problem as a 
constrained optimization problem 


(ero) — aus (5.22) 
The alternating direction method of multipliers (ADMM) (Glowinski 
and Marrocco, 1975; Gabay and Mercier, 1976) was introduced into the 
context of FFI-based methods by Michel et al. (2000; 2001), see also the 


application to non-smooth optimization by Willot (2020). Applied to our 
problem, we investigate the augmented Lagrangian function 


Dolg, ev) = F(E) + 9(e) + (v, £- e) +$ lE- ell’, (6.23) 


involving a penalty factor p > 0 and the Lagrange multiplier, i.e., our 
flow field, v : Yy — R®. The damped ADMM recursion (4.34) with 
damping factor ô is given by 


gk+t1/2 — argmin.L,(8, eo”), 
et = 2(1 _ dygktt/2 u (1 _ 25)e*, 
k+1 ; k+1 k (5.24) 
eT: = aremin, L(t" ev”), 


yktl = yk EB pier = gern), 
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Ux Ux 


Figure 5.2: Schematic of the projection of the vector v onto the admissible set in the case of 
isotropic (left) and anisotropic (right) crack resistance. (Ernesti and Schneider, 2022) 


More explicitly, the ADMM algorithm for anisotropic minimum cut/ 
maximum flow with adaptive penalty factor p” is given by 


gktl/2 = Pk 


EEH — 9(1 — E+ — (1 —28)e¥, (5.25) 
ektl = [v* + pe ger — Pee (v* + pr | /p*, 


yktl = yk ai pr ce Z gern). 


where Px, and Pea denote the orthogonal projectors onto the sets Kz 
and Ca, respectively. 


5.3.3 The anisotropic projector problem 


Within the ADMM iterations (5.25), evaluating two projection operators 
is required. The projection onto the compatible fields expressed via 
the projector Pk; may be efficiently performed with the help of the 
FFT (4.32). Additionally, the orthogonal projection onto the set Cg is 
required, expressed via the projection operator Pea and illustrated in 
Fig. 5.2. In the isotropic case, the set of constraints C comprises a sphere 
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per voxel. Thus, the projection onto C involves orthogonal projections 
onto a spheres with radii y(x) and is straightforward. In the anisotropic 
case, however, the set Ca comprises an ellipsoid for each voxel, with the 
eigenvalues y; of G as semi-axes. Following Kiseliov (1994), we write 
this projection as an optimization problem. 

Consider a vector v € R” and a symmetric, positive definite tensor Q. 
We seek the projection w € IR”, such that 


w-vl? + min . (5.26) 
wT Qw<1 


Introducing the Lagrange parameter X, the governing KKT-conditions, 
see for instance (Nocedal and Wright, 1999, Thm. 12.1), read 


2(w —v) +2AQw =0, (5.27) 
w Qu-1<0, (5.28) 
A>0 Aw? Qw-1)=0. (5.29) 


From conditions (5.27) we find 
w = (Id + AQ)7'v. (5.30) 


If the vector v satisfies v? Q v < 1, the problem (5.26) is trivially solved 
for w = v (and A = 0). We therefore focus on the case v? Q v > 0. Since 
the tensor Q is symmetric and positive definite, the inequality (5.28) 
describes a convex domain. Thus, the projection w lies on the boundary 
of the admissible set. Thus, the inequality (5.28) is satisfied as an equality. 
We insert the representation (5.30) into the conditions w? Qw — 1 = 0 
and solve the resulting equation 


uv? (Id +Q) Qld +Q) w -—1=0 (5.31) 
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for the scalar A, using Newton’s method and initial value \° = 0, follow- 
ing Kiseliov (1994). In the mentioned reference, global convergence of 
this algorithm is proved. 

For the application at hand, where n = 6, we set 


Q = G(x)? (5.32) 
for each x € Y. Therefore, the projection operator summarizes as 


v(2), eo) <1, 


Pegv(z) zu 
(Id+AQ)"!v(x), otherwise, 


with A solving (5.31). 


5.4 Numerical examples 


5.4.1 Setup 


The presented computational approach was integrated into an existing 
FFT-based code for computing effective crack energies on microstruc- 
tures (Ernesti and Schneider, 2021), which is embedded into an FFT- 
based computational homogenization code for thermal conductiv- 
ity (Dorn and Schneider, 2019). The software is written in Python with 
Cython extensions (and OpenMP). All computations were performed 
using the ADMM solver presented in chapter 4 with either the Barzilai- 
Borwein adaptivity or the averaging adaptivity, and a damping factor 
of 0.25. If not mentioned otherwise, the governing equations were 
solved for a prescribed tolerance tol = 1074 and the convergence 
criterion (Schneider, 2021b) 


Ls 2 
jet - ett? 


„Stel If): 6.33) 
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The computational experiments in this section were run on a desktop 
computer with 32GB RAM and six 3.7GHz cores and on a worksta- 
tion with 512 GB RAM and two Intel Xeon(R) Gold 6146 processors 
(12 x 3.20 GHz), respectively. 


5.4.2 A polycrystalline microstructure 


As our first example, we consider a polycrystalline microstructure gen- 
erated synthetically based on Laguerre tessellations and the algorithm 
described in Kuhn et al. (2020). Following a similar approach as Nguyen 
et al. (2017) within the context of phase-field fracture, we distinguish 
between 2D and 3D structures. In a 2D microstructure with distinct 
anisotropy, we identify one weak and one tough direction, whereas in 
3D, we identify one weak and two tough directions, resulting in one 
weak plane. Since elastic deformation and thus elastic material constants 
do not play a role in our model, we normalize the crack resistance tensor 
to a value y and consider a cleavage anisotropy factor 3 € [1,100] as 
proposed in Clayton and Knap (2014; 2015) for polycrystalline silicon 
carbide. The crack resistance tensor for 2D and 3D structures is given by 


G2 = ee | r A Rgrain and 

y o 0 (5.34) 
GP = Ran | 0 Py 0 | A, 

0 0 By 


respectively, where Rgrain is a rotation matrix encoding the orientation 
of each grain. 


Two-dimensional structures 


We start with a 2D Laguerre tessellation and planar grain orientations. 
In our first study, we vary the number of grains as well as the cleavage 
anisotropy factor 8. For each structure, we vary the loading angle 
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(0 p = 60° (d) ¢ = 90° 


Figure 5.3: Crack path through a microstructure with 256 grains and a cleavage anisotropy 
factor 8 = 10 for prescribed normal £ = (cos(y), sin()). (Ernesti and Schneider, 2022) 
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in seven equidistant steps from 0 to 90 degrees. For each loading 
angle, we additionally consider the case where each grain is rotated 
by 90 degrees. This results in 14 computations per structure. We used 
the Barzilai-Borwein adaptivity within the ADMM solver for most 
computations. For some cases, the averaging adaptivity converged 
faster and we switched to the latter solver choice in those cases. We 
extracted the mean value as well as the standard deviation of the 14 
computations per structure. The crack path for various loading angles 
is shown in Fig. 5.3. We see that the crack changes its direction for each 
grain in order to minimize its surface energy. To close the crack path for 
a non-axis aligned normal, the crack has to pass the cell several times. 
Furthermore, we observe local similarities in the crack path for different 
loading angles. 

Fig. 5.4a shows the influence of the number of grains on the effective 
crack energy. For a low number of grains, we observe a rather large 
standard deviation and no clear trend in the mean value. From 16? to 
32? grains, the standard deviation decreases significantly. This trend 
continues for 64” grains, where the standard deviation decreases even 
further while the mean value remains within the same range, indi- 
cating representativity. We thus find a structure of 32? grains to be 
sufficiently large. 

Secondly, we investigate the influence of the anisotropy factor on the ef- 
fective crack energy for the structure containing 1024 grains, see Fig. 5.4b. 
Note that the case 6 = 1 gives rise to a homogeneous £-field, since no 
local differences arise in the crack resistance. In particular, the effective 
crack energy becomes Yet = y. The effective crack energy increases with 
an increasing anisotropy factor until a certain threshold is reached at 
b = 50, beyond this internal contrast, no significant increase is visible. 
At this threshold, the crack favors the weakest direction in each grain. 
This is also reflected in Fig. 5.5. Clear differences in the crack path may 
be observed between 3 = 5 and $ = 10, see Fig. 5.5a and Fig. 5.5b, 
respectively. For increasing anisotropy, the differences become more 
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(a) Influence of the number of grains on the 
effective crack energy, 6 = 10 


(b) Influence of the anisotropy factor 8 on the 
effective crack energy for 1024 grains 


Figure 5.4: Influence of the number of grains and the cleavage factor on the effective crack 
energy in the two-dimensional case. (Ernesti and Schneider, 2022) 


subtle, such that the crack paths for 8 = 20 and 3 = 50 are almost 
indistinguishable, see Fig. 5.5c and see Fig. 5.5d, respectively. 


Three-dimensional structures 


We consider generated 3D Laguerre tessellations (Kuhn et al., 2020) 
with random grain orientation in each grain. Similar to the 2D case, 
we vary both the number of grains and the anisotropy factor. Since 
3D simulations are much more costly compared to 2D simulations, we 
perform only three simulations per microstructure for the size study by 
investigating a normals in ez, ey and e,-direction, respectively, and only 
one simulation per cleavage anisotropy factor 8 for the internal contrast 
study. Fig. 5.6 shows two cracked microstructures with 216 and 1728 
grains, respectively. Similar to the 2D case, the crack surface changes its 
orientation in each grain. Fig. 5.7a shows the results of the size study 
for an anisotropy factor 8 = 10. The effective crack energy increases 
for increasing number of grains. For a very small structure containing 
27 grains, the deviation between the three loading directions is rather 
large. For an increasing number of grains, this standard deviation is 
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(d) 8 = 50 


Figure 5.5: Crack path through a microstructure containing 1024 grains for € = ez and 
varying cleavage anisotropy factors. (Ernesti and Schneider, 2022) 
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(a) Cracked structure containing 216 grains (b) Cracked structure containing 1728 grains 


Figure 5.6: Cracked microstructure for different number of grains. (Ernesti and Schneider, 
2022) 


reduced to a negligible amount. Two main distinctions between the 
2D and the 3D case emerge in case of the size study. Firstly, the range 
of the effective crack resistances differ significantly: In the 2D case the 
effective crack energy ranged around 7° ~ 1.47, whereas for the 3D 
case we find 7°! > 5.2 y for 8 = 10. Secondly, the effective crack energy 
increases with increasing number of grains in the 3D, at least within the 
range of our observation, whereas in the 2D case we observed saturation. 
This indicates that even more than the considered 13824 grains could be 
necessary to reach representative results. 

Shifting our focus to the contrast study, see Fig. 5.7b, we observe a nearly 
linear correlation between the cleavage factor 3 and the effective crack 
energy. This clearly differs from the 2D case, which displayed a satura- 
tion. This effect has a geometric origin. In the 2D case, a continuous crack 
path can be found which passes each grain in the energetically most 
favorable direction. In 3D, this is not the case, as planar cracks within 
each grain cannot be combined into a continuous, global crack surface, in 
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(a) Influence of the number of grains on the (b) Influence of the material contrast on the 
effective crack energy, 8 = 10 effective crack energy for 1024 grains 


Figure 5.7: Influence of the number of grains and of the material contrast on the effective 
crack energy in the three-dimensional case. (Ernesti and Schneider, 2022) 


general. Hence, in the two- and three-dimensional context, the effective 
crack resistance of polycrystalline materials differs fundamentally in the 
context of anisotropic intergranular fracture, modeled with one cleavage 
parameter. In two spatial dimensions, the parameter 5 plays a subordi- 
nate role (at least if it is sufficiently large), in the three-dimensional case 
we see clearly that 3 is an important material parameter, which needs to 
be identified. 


5.4.3 A carbon fiber reinforced polymer 


In this section, we investigate a carbon fiber reinforced composite. Car- 
bon fibers are used due to their favorable strength-density ratio. The 
individual fibers have a strong anisotropy in both elastic and strength 
properties. As a result of the manufacturing process, they show higher 
Young’s modulus and higher strength in fiber direction compared to the 
plane perpendicular to it. We use a crack resistance tensor for a fiber 
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(c) Crack surface for € = ey, Y° = 0.85 y (d) Crack surface for € = ez, 7° = 0.877 


Figure 5.8: Carbon fiber reinforced composite and crack surface for varying direction. 
(Ernesti and Schneider, 2022) 
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oriented in (unit) direction p as 
G = 25y7p@p+0.547(Id—p 8 p), 


assuming an internal contrast of 50 for the crack resistance, as suggested 
by Pillai et al. (2020). Furthermore, we model the matrix material in 
our composite with the isotropic crack resistance y, assuming that the 
fibers transverse crack resistance is lower than that of the matrix material. 
The microstructure under consideration contains 290 fibers of 450um 
length and 7m diameter, oriented in an almost unidirectional manner in 
x-direction with second order orientation tensor diag(0.9, 0.05, 0.05) and 
a total of filler content of 15%, see Fig. 5.8a. The structure was generated 
using sequential addition and migration (Schneider, 2017). 

Fig. 5.8b shows the crack surface oriented in x-direction. We notice fiber 
pullout and matrix failure, as well as fiber damage that appears to be 
non-perpendicular to the fiber direction. The effective crack energy is 
increased by 50% compared to the matrix material. The crack surfaces in 
y- and z-direction are shown in Fig. 5.8c and Fig. 5.8d, respectively. Both 
crack surfaces look qualitatively similar. We notice both matrix failure 
and inter-fiber debonding, since we assumed the fibers perpendicular to 
their orientation to be weaker than the matrix material. For both cases, 
the effective properties are lower compared to the matrix material and 
almost equal. 


5.5 Conclusion 


In this chapter, we presented a numerical method for computing 
the effective crack energy of a heterogeneous medium with distinct 
anisotropy via weighted minimal surfaces. We derived the anisotropic 
maximum flow problem and discussed the implementation into an 
FFT-based homogenization tool for isotropic fracture. We saw that 
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both the solver framework and the discretization established for the 
isotropic case serves as a firm foundation for the anisotropic case. Indeed, 
the extension requires evaluating the projection onto ellipsoids in an 
efficient manner. 

This anisotropic extension of the homogenization of the fracture energy 
enables the treatment of additional material classes compared to 
previous works. Applications were presented for polycrystalline 
ceramics and carbon-fiber reinforced composites. 

In the literature on anisotropic phase-field fracture, 2D polycrystalline 
materials are often investigated in addition to the 3D case. The 
anisotropy may be encoded by the cleavage parameter 6, which 
penalizes crack propagation in certain directions and forces the crack 
path to follow a weak plane. In our homogenization framework, we 
observed the behavior for the 2D case to fundamentally differ from the 
3D case. In two dimensions, a crack is not geometrically restricted to 
follow the weakest plane through each grain. Therefore, the cleavage 
parameter 3 has no further meaning beyond a certain threshold. In the 
3D case, on the other hand, neighborhood relations between different 
grains prohibit the crack to form freely in order to follow the weakest 
plane for each grain. Therefore, we observe a strong dependence of 
the effective crack energy on the cleavage parameter 3, emphasizing its 
importance from the viewpoint of materials science. 

Additionally, we saw that our framework allows us to investigate carbon 
fiber structures, which show distinct anisotropy within each fiber. This 
enables modeling additional effects compared to isotropic fibers, studied 
in chapter 4. With isotropic inclusions either weakening, similar to pores, 
or toughening with respect to the matrix material can be modeled. The 
inclusion of anisotropic fibers allows for toughening in one direction, for 
instance the fiber direction, and weakening in the transverse direction, 
broadening the spectrum for material design. 
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Investigations on the influence of 
the boundary conditions when 
computing the effective crack 
energy of random heterogeneous 
materials using fast marching 
methods’ 


6.1 Introduction 


Nearly 25 years ago, Francfort and Marigo introduced a variational 
fracture model (Francfort and Marigo, 1998) based on Griffth’s original 
energy criterion (Griffith, 1921). The model was formulated in terms of 
an energy functional to be minimized and added impetus to phase-field 
fracture mechanics (Bourdin et al., 2000; Wu et al., 2020). The latter is 
rather popular, as it permits to simulate both crack initiation and 
crack propagation in a common framework based on traditional 


1 This chapter is based on Ernesti et al. (2023) which was built on the results of the master 
thesis by Lendvai (2022) which was supervised by me. In order to include this paper 
into the structure of this work I shortened both the introduction and the second section. 
Furthermore, minor additional changes to the manuscript have been made. The results 
which were originally presented in Lendvai (2022) are found in section 6.4. Figures 
which are found in the master thesis are labeled accordingly. 
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finite element methods. The phase-field model may be interpreted 
as a regularization of the Francfort-Marigo model, inspired by the 
Ambrosio-Tortorelli approximation (Ambrosio and Tortorelli, 1990) of 
the Mumford-Shah functional (Mumford and Shah, 1989). Alternatively, 
the phase-field fracture model may be interpreted as a non-local damage 
model (Dimitrijevic and Hackl, 2008; Bazant, 1991). 

Multi-scale methods (Milton, 2002) are used to predict the macroscopic 
behavior of microstructured materials, explicitly accounting for the 
material behavior of the constituents and their arrangement on the 
microscale. Homogenization theories serve as the mathematical un- 
derpinnings of modern multi-scale methods. Seminal contributions 
focus on periodic homogenization (De Giorgi and Spagnolo, 1973; 
Babuska, 1973; Larsen, 1975), where the underlying microstructure 
is given by a periodic cell. Since microstructured materials are often 
random due to their manufacturing process (Jeulin, 2021), stochastic 
homogenization results (Papanicolaou and Varadhan, 1981; Kozlov, 
1978) have been established, where an infinitely large domain with a 
random microstructure is investigated. One method to evaluate effective 
properties based on stochastic homogenization results is by exploring 
representative volume elements (RVEs). For prescribed accuracy and fixed 
boundary condition type, a volume element is called representative 
provided it approximates the effective properties of the infinitely large 
domain up to this prescribed accuracy (Drugan and Willis, 1996). Since 
the size of the RVE is not known a priori and dependent on various 
factors, the RVE size is often found via computations on volume elements 
of increasing size (Gusev, 1997; Segurado and Llorca, 2002) until a 
desired degree of representativity is reached. For elliptic PDEs such as 
elasticity and conductivity problems, theoretical results show that for 
different boundary conditions applied to cells of finite size, i.e., periodic, 
Dirichlet or Neumann boundary conditions, the effective quantities 
evaluated on the cells share the same infinite-volume limit (Sab, 1992; 
Bourgeat and Piatnitski, 2004; Owhadi, 2003). Yet, it is well known that 
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the boundary conditions may have a strong influence on the necessary 
size of the RVE, see Kanit et al. (2003). 

In the field of fracture mechanics, multi-scale approaches face additional 
difficulties compared to linear elastic or conductivity problems. One 
necessary ingredient for multi-scale methods is a distinct scale sepa- 
ration, allowing the quantities of interest, displacements, stresses or 
strains, to be separable into a large-scale and small-scale component. In 
an elastic material without a crack, for instance, this separation leads to 
a cell formula on the microscale, from which effective quantities may be 
derived. This approach does not work within a cracked microstructure, 
since this crack and the stress singularity resulting from it would be 
present on both the micro- and the macroscale. On the other hand, 
investigations on volume elements of finite size may be conducted 
with phase-field fracture models (Chen et al., 2019; Ernesti et al., 2020). 
However, it is well known (Gitman et al., 2007) that such a procedure 
will not, in general, lead to a macroscopic model for softening materials. 
Thus, special care is required. 

Our strategy is based on a periodic homogenization result of Braides et al. 
(1996) for the Mumford-Shah functional (Mumford and Shah, 1989) and 
energetic minimization. This result covers the Francfort-Marigo model 
of brittle fracture (Francfort and Marigo, 1998) in the case of anti-plane 
shear when considering a fixed quasi-static time discretization and 
neglecting crack irreversibility (for instance in the very first load step 
on an un-cracked specimen). Within their result, Braides et al. (1996) 
show a decoupling of the volumetric part and the surface part of the 
Mumford-Shah functional. In the context of the Francfort-Marigo model 
this implies a decoupling of the effective stiffness and the effective crack 
energy in the anti-plane shear case. Furthermore, they provide specific 
formulas for both effective quantities. From their work the effective crack 
energy is defined as the area of the crack resistance-weighted minimal 
surface cutting through the microstructure. Numerical approaches to 
computing the effective crack energy have been proposed by Schneider 
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(2020), as well as Ernesti and Schneider (2021; 2022) using FFT-based 
algorithms and periodic boundary conditions, see also chapters 4 and 5. 
Recently, Michel-Suquet addressed the approach based on the ho- 
mogenization result of Braides et al. (1996) using an alternative so- 
lution strategy by pointing out similarities with limit load analy- 
sis (Christiansen, 1981). 


Contributions 


The aim of this chapter is to investigate the effective crack energy of 
solids with random microstructures and the influence of the imposed 
boundary conditions. We give a brief summary on the cell formula for 
the effective crack energy based on the periodic (Braides et al., 1996) and 
stochastic (Cagnetti et al., 2019) homogenization results for the Mumford- 
Shah functional (Mumford and Shah, 1989) in section 6.2. For both the 
periodic and the stochastic setting, the effective crack energy is expressed 
in terms of a multi-cell formula on specifically notched cubes, which 
we call Dirichlet boundary conditions. For periodic homogenization, 
this multi-cell strategy is overly arduous, and a single-cell formula 
is sufficient, provided periodic boundary conditions are used, as well. 
Therefore, it makes sense to investigate periodic boundary conditions 
for random materials, as well. Indeed, for periodic boundary condi- 
tions, powerful numerical tools for computing the effective quantities 
based on the fast Fourier transform (FFT) are available (Schneider, 2020; 
Ernesti and Schneider, 2021; 2022), see also chapters 4 and 5. Dirichlet 
boundary conditions on the other hand are not easily integrated into 
this framework. 

In a two-dimensional setting the problem of computing the effective 
crack energy reduces to finding weighted minimal paths. One prominent 
algorithm to compute such paths is the fast marching method intro- 
duced by Sethian (1996; 1999) where fast implementations are publicly 
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available’. In the context of fracture mechanics, the fast marching 
method has already been used for fatigue fracture using stress intensity 
factors (Jovičić et al., 2005) and in combination with the extended finite 
element method (Stolarska et al., 2001; Sukumar et al., 2003; 2008). We 
discuss a straightforward way to compute the effective crack energy with 
the help of the fast marching method in section 6.3. One advantage of 
this method is that Dirichlet boundary conditions can easily be applied, 
since every point of a domain may be used as a starting or ending point. 
Section 6.4 comprises the numerical results. We first validate the fast 
marching method in terms of accuracy of the discretization and compare 
the results for periodic boundary conditions with established tools. 
Finally, we investigate the influence of the boundary condition on the ap- 
proximated effective crack energy for microstructure samples of increas- 
ing cell size. We compare Dirichlet and periodic boundary conditions 
and study their necessary size of the computational cell. 


6.2 The effective crack energy of heteroge- 
neous random media’ 


Consider a domain Q C R and a heterogeneous field of crack resistances 
y : Q — Rso which is periodic with periodicity 7. Following the 
homogenization result of Braides et al. (1996) for the Mumford-Shah 
functional (Mumford and Shah, 1989) we define the effective crack 
energy as follows: Inside an infinite periodic continuation of our material 
with periodicity 7 we place a cube LQ,, of edge length L with its eı axis 
rotated onto the prescribed normal n. On such a cube, we compute a 
y-weighted minimal surface $ cutting the cube under the constraint 


2 https:/ /github.com/scikit-fmm/scikit-fmm, accessed in November 2021 
3 In this section we provide a brief summary of section 2.4.2 to enhance the readability of 
this chapter. 
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that the surface cuts the boundary of the cube at x; = L/2 within the 
coordinate system of the cube, see Fig. 2.10 in section 2.4.2 for a visual 
representation. The effective crack energy is given by the limit of these 
computed weighted minimal surfaces as the cube edge-length L — oo. 
In mathematical terms, the effective crack energy is defined via 


m inf 


. 1 
= im n gaa [144 (6.1) 


Yett(7) 


Let us take a closer look at the cell formula (6.1). From a computational 
point of view, the limit of L — oo is not practicable, as we can only deal 
with finite computational domains. To overcome this issue, we may 
restrict to a single cell LQ and employ the boundary conditions used 
above, which we call Dirichlet boundary conditions. In this case, we fix 
the surface S on the boundary OLQ of the cube at xı = L/2 within the 
local coordinate system of the cube. Actually, any cut at x; € [0, L] could 
be chosen. We choose x; = L/2 for definiteness. 

The approach by Braides et al. (1996) involves a multi-cell formula 
(6.1) although the homogenization problem is periodic. As the inte- 
grand in the problem (6.1) is convex, it is reasonable to hope that a 
single-cell formula may prove sufficient for computing the effective 
crack resistance er. This is indeed true, as shown by Braides and 
Piat (1995) and Chambolle and Thouroude (2009). More precisely, they 
showed that the effective crack energy Yer(n) in equation (6.1) may be 
computed on a single cell Y, with period 7 and for fields with periodic 
boundary conditions 


1 
eH) = oy, Peace, Met Valle 62 


Many materials of industrial relevance show a randomness in their mi- 
crostructure composition, and periodic homogenization is not sufficient 
for describing those. Recently, Cagnetti et al. (2019) proved an extension 
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of the result of Braides et al. (1996) to stochastic homogenization for 
the Mumford-Shah functional. Remarkably, the result provides a true 
extension of the periodic case. The explicit formula they find is the 
same as in the periodic case, i.e., (6.1), but the infinite-volume limit is 
indispensable for stochastic homogenization. 

The case of stochastic homogenization of elastic materials is well-studied. 
To compute effective quantities on cells of finite size, appropriate bound- 
ary conditions are required. Upon the infinite-volume limit the chosen 
boundary conditions do not affect the effective quantities (Sab, 1992; 
Bourgeat and Piatnitski, 2004; Owhadi, 2003). However, for cells of finite 
size, the boundary conditions have an influence on the approximation 
quality of the "true" effective stiffness, as shown by Sab (1992); Kanit et al. 
(2003). Typically, using periodic boundary conditions and periodized 
ensembles of random microstructures result in optimal convergence 
rates (Schneider et al., 2022). 

In contrast to computing the effective stiffness or effective conductivity 
of a microstructured material, very little is known about the influence 
of the boundary conditions when computing the effective crack energy 
evaluated on cells of finite size. The aim of this chapter is to provide a 
first step in this direction. For the periodic boundary conditions and in 
three spatial dimensions, Schneider (2020) proposed an algorithm for 
computing the effective crack energy on cells on finite size. The approach 
relies on Strang’s minimum cut/maximum flow duality (Strang, 1983), 
see also chapters 4 and 5. 

Due to the tremendous computational effort involved, we restrict to 
two-dimensional microstructures. In this case, it is possible to compute 
shortest paths with fast marching, see section 6.3, which is well-known 
among experts. 
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6.3 Finding the minimum cut with the fast 
marching method 


6.3.1 From minimum cut to shortest path problems 


J 


r 


(a) Microstructure (b) ||€|| for € = e1 (c) ||€|| for 
E = (3e1 + e2)/ V10 


Figure 6.1: Microstructure and minimum cut field (€ in eq. (2.32)) for an axis-aligned and 
a non-axis-aligned prescribed mean normal. (Ernesti et al., 2023) 


To gain some intuition into minimum cut fields in two spatial dimen- 
sions, we consider a periodic microstructure of circular inclusions, 
shown in Fig. 6.la. The structure contains 35 fillers, i.e., 50% area 
fraction. We consider tough inclusions, i.e., these have a (much) higher 
crack resistance than the matrix material. Fig. 6.1b and Fig. 6.1c show 
the minimum cut for mean crack normals € = e; and € = (3e1 + e2)/ VTO. 
In both cases, the minimum cut field localizes and takes the shape of 
a crack path that cuts through the microstructure. However, a distinct 
difference is present between the two cases. For the axis-aligned case, the 
cut field traverses the microstructure once, cutting from left to right. For 
the non-axis aligned case, the cut field wraps around the microstructure 
several times in order to both preserve the mean normal of the cut and 
to retain periodicity. 

Thus, at least for the axis-aligned case, it appears reasonable the 
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minimum cut may be computed by an algorithm which returns a 
(weighted) shortest path (Jeulin, 1988; 1994a; Noyel et al., 2011). In- 
deed, after fixing two corresponding points on opposing edges of the 
microstructure, a minimum weighted path joining the two points would 
have to be computed. 


QL 
(a) Obstacles with loop holes (b) Diagonal obstacle 


Figure 6.2: Periodic microstructures LQ (with a periodic extension in y-direction), giving 
rise to possible problems for shortest path methods. (Ernesti et al., 2023) 


Some caution is advised with this kind of strategy, as the following 
two examples demonstrate. For a start, we consider the periodic 
microstructure shown in Fig. 6.2a with imposed crack normal & = cı. 
If the crack resistance in the rectangles (shown in blue) is much higher 
than in the complement (shown in white), the minimum cut is forced 
to navigate through the white pathways. In this process, more than 
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one unit cell needs to be crossed. In the example shown, the green 
curve crosses the horizontal "boundary" twice. Such a curve may be 
represented by a shortest path algorithm if periodicity in y-direction is 
accounted for. Otherwise the red curve would arise as the shortest path 
from left to right. 

Unfortunately, taking periodicity in y-direction into account does not 
always offer the proper strategy, as the microstructure in Fig. 6.2b 
shows. A straight "obstacle" with high crack resistance is placed along 
the diagonal. For prescribed normal E = e, the minimum cut has to 
cross the obstacle. The shortest path strategy with periodic boundary 
conditions in y-direction, however, would give rise to the green path. 
Unfortunately, the shown path does not give rise to the correct path 
normal E = eı, but to the normal n pointing in diagonal direction! If, 
instead, no periodicity in y-direction is permitted, the correct crack path 
(in red) is computed. 

To summarize, the charming idea of working with shortest path 
algorithms to compute the effective crack energy for axis-aligned crack 
normals may be unsuited to some microstructures. Therefore, it is 
unavoidable to perform a validation against minimum cut methods. For 
the microstructure models considered in this article, such a comparison 
is contained in section 6.4.3. 


6.3.2 Finding shortest paths by the fast marching method 


There is a deep connection between the eikonal equation and efficient 
path finding which led Sethian (1996; 1999) to devise computationally 
efficient algorithms for the latter. More precisely, consider a domain 
Q c R“ and suppose a wave propagates through our domain, starting 
from some point zo € Q at a given velocity v : Q — Ryo. Then, the 
time T : Q — Rso this wave needs to arrive at point x € Q solves the 
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eikonal equation 
IVT) =; we, (6.3) 


with T (xo) = 0. If the velocity v is spatially homogeneous, the level sets 
of the travel time T describe concentric spheres around the starting point 
xo. For a heterogeneous velocity, the wavefront is refracted. Sethian 
(1996; 1999) introduced the fast marching method as a fast algorithm for 
solving the eikonal equation. The fast marching method is an integrated 
strategy where the spatial discretization and the strategy for solving the 
eikonal equation are well orchestrated. More precisely, the solution 
strategy uses a modification of Dijkstra’s algorithm (Dijkstra, 1959) 
well-known in graph theory (Gera et al., 2018). The fast marching 
method finds application in various fields ranging from shortest path 
finding (Kimmel and Sethian, 1996; Mirebeau, 2018; Garrido et al., 2006) 
to simulating wildfire spreading (Carballeira et al., 2021) and within the 
extended finite element method (Sukumar et al., 2003; 2008). 

In d spatial dimensions and on a regular grid with N d grid points, the 
fast marching method has the computational complexity O(N“ log N). 
In contrast, iterative procedures for solving the eikonal equation (6.3) 
typically have a complexity of O(N“*!). This complexity reduction is 
partly caused by an underlying min-heap data structure (Williams, 1964). 
The problem of computing the effective crack energy on a microstructure 
involves finding a weighted minimal surface, as we pointed out in 
chapter 2. For the special case of two-dimensional structures, this 
problem simplifies to finding shortest paths in a given two-dimensional 
microstructure, for which various methods are available Willot (2015b; 
2019). In particular, the fast marching method may be applied as follows. 


1. The crack resistance y(x) serves as the weight in computing the 
weighted shortest path playing the role of a resistance for the crack to 
propagate. In contrast, for a propagating wave, the velocity v(x) 
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£ 


Figure 6.3: Example for (non-unique) shortests path with Dirichlet (index "D", red) and 
periodic (index "P", green) boundary conditions. (Ernesti et al., 2023) 


enables the propagating wave to travel faster at a higher speed. 
Therefore we set the right hand side of the eikonal equation (6.3) 
to 7(x) instead of 1/v(x) for computing the effective crack energy 
with fast marching". 

2. The solution field T(x) embodies the y-weighted distance from point 
x to the origin x°. We therefore call it the distance field throughout 
this work. 

3. To compute the effective crack energy of a given microstructure 
cell Y = [0, L]? with crack normal n = e and Dirichlet boundary 
conditions, we choose a starting point x = (0, L/2)T and evaluate 
the distance field T in x* = (L, L/2)”. The effective crack energy is 
given by T(x*)/L, see the red path in Fig. 6.3. 


4 The crack resistance y and the inverse velocity 1/v have different physical units. 
However, both the effective crack energy and the travel time field scale homogeneously 
under a rescaling of the crack resistance and the inverse velocity. Thus, upon introducing 
a conversion factor between crack resistance and the velocity in the beginning, the same 
conversion factor permits to recover the effective crack energy from the computed 
distance field. For simplicity of notation, we therefore suppress mentioning the 
conversion factor and tacitly assume it to be chosen appropriately. 
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LL. YR 


x 


(a) Dirichlet bc (b) Mixed bc (c) Periodic bc 


x 


Figure 6.4: Distance field and crack path for different boundary conditions for a single 
circular inclusion microstructure. (Ernesti et al., 2023) 


4. The solution for periodic boundary conditions with mean normal 
n = e, is given by the y-weighted shortest periodic path from the 
left hand side to the right hand side. To find this shortest path with 
the fast marching approach, we consider all paths with the starting 
point x° = (0,y) and end point x* = (L,y) for some y € [0, L], and 
select the starting point which returns the smallest crack energy, see 
the green path in Fig. 6.3. On a computational grid with N x N pixels, 
this process includes N fast marching computations, which increases 
the computational complexity to O(N? log N). 

5. The fast marching method enables computing the minimum crack 
energy in a straightforward way. However, the involved crack path is 
not directly accessible. Rather, different approaches are available to 
obtain the crack in post processing, see for instance Noyel et al. (2011). 
Using the fact that the shortest crack path is perpendicular to the level 
sets of the distance field T, we rely on a gradient descent method 
to compute the crack path from any point z to the origin x. To do 
so, we compute a spline interpolation on the numerically evaluated 
gradient of the distance field T. 

One advantage of the fast marching method over the maximum flow 

approach discussed in chapters 4 and 5 is that additional boundary 
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conditions can easily be studied, as we can choose any point of Q as our 
starting or ending point. We consider three different cases, illustrated 
in Fig. 6.4, which shows the level sets of the distance field and the 
resulting crack path for a structure containing a single circular inclusion 
of diameter L/2 positioned at the center of a square with edge length L. 
The inclusion has a much higher crack resistance than the embedding 
matrix, forcing the crack path to avoid the inclusion altogether. The 
distance field and the crack path for Dirichlet boundary conditions is 
shown in Fig. 6.4a. The distance field describes concentric circles starting 
from the left hand side until the inclusion is reached and a refraction 
occurs. To draw the crack path, we start on the right hand side at y = L/2 
and follow the path perpendicular to the contour lines of the distance 
field. Notice that this path does not prescribe the shortest path from the 
right hand side of the structure to the origin on the left hand side. Indeed, 
Fig. 6.4b shows this shortest path originating in x° to the right-hand side. 
To draw this path we select x* as the point on the right hand side with 
the minimum effective crack energy and draw the path perpendicular to 
the contour lines of the distance field. Since the mean crack normal may 
differ from the prescribed mean normal we do not focus on this case 
in this paper. For periodic boundary conditions the crack path is not 
unique for this example, as both above and below the inclusion straight 
paths are possible, see Fig. 6.4c for one possible option. 
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6.4 Numerical investigations 


6.4.1 Setup’ 


The fast marching based algorithm for computing the effective crack en- 
ergy was implemented in Python 3 based on the scikit-£mm module®, 
which provides both first-order and second-order fast marching methods. 
These differ in the convergence order of the underlying approximation 
of first derivatives. The first-order fast marching method uses classical 
forward and backward differences, whereas for the second-order method 
a three point stencil approximation of forward and backward differences 
is used Rickett and Fomel (1999). The crack paths were visualized by 
first computing the gradient of the resulting distance field by finite 
differences, interpolating this gradient field with bi-cubic splines and 
finally using gradient descent starting from the end point of the crack 
path. 

For the computations based on the minimum cut/maximum flow 
approach, we relied on an in-house FFT-based code (Schneider, 2020; 
Ernesti and Schneider, 2021; 2022), see also chapters 4 and 5. The contin- 
uous equations were discretized with the CCMF discretization (Couprie 
et al., 2011) and solved by a damped version of the alternating direction 
method of multipliers (ADMM) (Glowinski and Marrocco, 1975; Gabay 
and Mercier, 1976) with adaptive choice for the penalty parameter, see 
Ernesti and Schneider (2021). We chose a relative tolerance of 1074. 

All fast marching computations were run on an ARM-based SoC 
Apple M1 with 8 GB of RAM using a single thread. The minimum 
cut/maximum flow computations were performed on a desktop com- 


5 The numerical investigations presented in this section are the result of the master thesis 
by Lendvai (2022) which was supervised by me. We reevaluated the findings of this 
thesis in a joint publication (Ernesti and Schneider, 2022) which forms the basis of this 
chapter. 

6 https: //github.com/scikit-fmm/scikit-fmm, accessed in November 2021 
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y L R 


(a) Microstructure and (b) Distance field for Dirichlet boundary (c) Resulting crack path 
anticipated crack path conditions 


Figure 6.5: Microstructure of a rotated square with anticipated crack path, distance field 
and crack path for Dirichlet boundary conditions. (Lendvai, 2022; Ernesti et al., 2023) 


puter with 32 GB of RAM and six 3.7 GHz cores. 


6.4.2 A single rotated square inclusion 


To investigate the accuracy of the fast-marching approach, we start with a 
structure containing a single square inclusion of edge length L/2, which 
is rotated at 45 degrees and positioned at the center of our computational 
domain of edge length L, see Fig. 6.5a. The crack resistance of the 
inclusion is given by 7b = 10 Ymat- We consider an initial double notch 
crack at y = 0.5 L, i.e., Dirichlet boundary conditions. The analytical 
solution, which may be extracted from the structural measurements, see 
Fig. 6.5a, is Yett/Ymat = Y3/2 ~ 1.22247, see Fig. 6.5a. Fig. 6.5b shows 
a contour plot of the distance field. Starting from the initial notch on 
the left hand side, the distance field initially shows a circular expanding 
front. Upon hitting the inclusion, the field is refracted and a new circular 
front with the top/bottom edge of the inclusion as origin continues to 
the other side. Since the distance field is symmetric with respect to the 
x-axis, we slightly perturb the starting point for our gradient descent 
method in y-direction to break this symmetry and enforce a unique crack 
path. The evaluated crack path is shown in Fig. 6.5c. We observe that it 
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(a) Effective crack energy vs. number of pixels (b) Relative error of the effective crack energy vs 
per dimension number of pixels per dimension 


Figure 6.6: Effective crack energy and relative error for the rotated square microstructure. 
(Ernesti et al., 2023) 


matches the geometrically anticipated path. 

Next, we investigate the quality of the solution with respect to the grid 
size. We compute the effective crack energy for different resolutions 
ranging from 100 to 6400 pixels per side length of our computational 
domain. Furthermore, we investigate the performance of both first-order 
and second-order fast marching methods. The results are shown in 
Fig. 6.6a, where the absolute values of the crack energy are depicted, 
as well as in Fig. 6.6b, which shows the relative error compared to 
the analytical solution. Both the first and the second-order methods 
converge to the analytical solution with a linear rate of convergence. 
However, the second-order approach leads to a higher accuracy than the 
first-order approach, even on a coarser grid. To reach the accuracy of the 
second-order fast marching using first-order requires almost four times 
more pixels per edge length. Hence, we rely on the second-order fast 
marching method for the remainder of this chapter. 
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6.4.3 Comparison with minimum cut/maximum flow 


In our next study, we investigate a periodic structure containing 32 
unidirectional fibers and a filler fraction of 50%, which are represented 
as circular inclusions, see Fig 6.7. The structure was generated using 
a mechanical contraction algorithm (Williams and Philipse, 2003). The 
inclusions are considered tough, i.e., they have a higher crack resistance 
than the matrix and we set yab = 10 Ymat- We investigate various 
resolutions ranging from 128 to 4096 pixels per edge length. 

In this section, we would like to investigate whether the results obtained 
with the fast-marching method and with the FFT-based technique give 
rise to similar results. Indeed, as discussed in section 6.3.1, it is necessary 
to exclude certain pathological situations in order to gain confidence 
into the results obtained with the fast marching method. We wish to 
emphasize that due to the lack of uniqueness of the solutions to the 
minimization problem (2.32), we may only expect the obtained effective 
crack energies to be close. However, a similar obtained minimum cracks 
are certainly a sufficient condition for the latter. 

To enforce periodic boundary conditions in the fast-marching setting, 
we iterate over all pixels on the right hand side of the microstructure, 
evaluate the distance field at the same height on the other side and select 
the minimum. The distance field is shown in Fig. 6.7c. From the near 
center of the y-axis we observe a circular expanding front which is re- 
fracted at every inclusion. For both, the fast marching and the minimum 
cut/maximum flow formulation and the two considered resolutions, the 
crack paths are shown in Fig. 6.8. Notice that the way these two methods 
extract the crack path is very different. Whereas the crack path of the fast 
marching method is computed via gradient descent along the distance 
field, the crack path of the minimum cut/maximum flow approach is 
given by the total minimum cut through the microstructure with mean 
normal e,, which is a field that localizes around the crack attaining large 
values whose magnitude has no physical meaning. This results from 
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0.0 
y (a) Microstructure (b) Microstructure (c) Distance field for 
* discretized ona 128° grid discretized ona 512? grid periodic boundary 
conditions 


Figure 6.7: Microstructure containing 32 circular inclusions for different resolutions and 
distance field for periodic boundary conditions. (Lendvai, 2022; Ernesti et al., 2023) 


the fact that the minimum cut is given by the gradient of the periodic 
field ¢ in equation(6.2) which has a jump discontinuity across the crack. 
Evaluating this quantity numerically results in large but finite values 
which tend to infinity as the pixel length goes to zero. Both methods 
find extremely similar crack paths. On a coarser grid of 128? pixels, the 
fast marching crack in Fig. 6.8b exhibits some small isolation distance 
to the inclusions, which is not the case for the minimum cut field. This 
isolation distance vanishes for higher resolutions, see Fig. 6.8d. These 
computational results resolve possible doubts about the expressivity of 
the fast marching results for the considered microstructures. 

In Fig. 6.9, the effective crack energy for the two approaches under 
consideration is plotted against the resolution, together with the relative 
error, where we used the solution on the 4096? grid as the ground truth 
for each method. Both approaches show a linear rate of convergence with 
respect to the resolution per edge length. Furthermore, both approaches 
overestimate the effective crack energy on a coarser grid. However, 
in order to reach the accuracy of the minimum cut, the fast marching 
method requires between 1.5 and 2 times the resolution of the edge 
length. Notice the difference in the complexity of the two methods. The 
complexity of the minimum cut/maximum flow is mainly driven by the 
FFT, which has a complexity of O(N? log N), as well as the number of 
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(© Minimum cut on a 512? grid (d) Fast marching crack path on a 512? grid 


Figure 6.8: Periodic crack paths for minimum cut/maximum flow and fast marching 
method of a microstructure containing 32 circular inclusions for different grid sizes. 
(Lendvai, 2022; Ernesti et al., 2023) 
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Figure 6.9: Effective crack energy and relative error comparing fast marching with 
minimum cut/maximum flow. (Ernesti et al., 2023) 


required iterations, which range between 2000 and 7000 in order to reach 
the desired accuracy of 1074. The fast marching approach has a com- 
plexity of O(N? log N) for periodic boundary conditions. As a result we 
noticed for resolutions below N = 2048 that the fast marching method 
required less computational time than the minimum cut/maximum 
flow method, each running on a single thread. For N = 4096 the fast 
marching method required twice the computational time of minimum 


cut/maximum flow. 


6.4.4 The influence of boundary conditions 


In our next study we investigate the influence of the boundary condi- 
tions on the effective crack energy for computational cells of increasing 
size. We consider two types of boundary conditions, namely Dirichlet 
boundary conditions and periodic boundary conditions. For Dirichlet 
boundary conditions, we consider the crack propagating at y = 0.5 L to 
the other side of the domain at the same height. Fully periodic boundary 
conditions are attained by the minimum value when iterating over 
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Figure 6.10: Crack paths for different boundary conditions. (Lendvai, 2022; Ernesti et al., 
2023) 


all pixels in y-direction using Dirichlet boundary conditions starting 
from each pixel. The material parameters are chosen as before, i.e., the 
inclusions are considered tough with a material contrast of 10. 

A comparison of the crack paths for different boundary conditions 
is shown in Fig. 6.10, where we consider a microstructure with 50% 
circular inclusions, i.e., unidirectional continuous fibers. The crack 
path for the periodic boundary conditions interacts with less inclusions, 
resulting in a path with more straight segments compared to the Dirichlet 
boundary conditions. 

To further investigate the boundary conditions, we consider microstruc- 
tures with 30% and 50% filler fraction and a varying number of inclu- 
sions ranging from 5? to 80°. For each number of inclusions we consider 
100 microstructure realizations which were generated using mechanical 
contraction (Williams and Philipse, 2003). For the Dirichlet boundary 
conditions we consider all realizations. To reduce the computational 
costs for periodic boundary conditions we only take half of the real- 
izations into account for a fiber count of 50° and higher. The results 
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Figure 6.11: Comparison of the boundary conditions for 30% filler content. (Ernesti et al., 
2023) 
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for volume fractions 30% and 50% are shown in Fig. 6.11 and Fig. 6.12, 
respectively. Fig. 6.1la and Fig. 6.12a show a histogram of the crack 
energy for 25? inclusions. On the y-axis, the number (in percent) of 
microstructures is shown whose effective crack energy corresponds to 
the x coordinate. We notice that the range of the periodic boundary 
condition is shifted to the lower values of the effective crack energy 
ranging, into the lower part of the Dirichlet boundary conditions. For the 
Dirichlet boundary conditions we notice some accumulation in the lower 
range up to Yes = 1.025 Ymat for 30% filler fraction and Yer = 1.06 Ymat 
for 50% filler fraction. Above these thresholds both histograms show 
some dispersion. These dispersions result from the fact that for some 
microstructures, the initial crack in the Dirichlet boundary conditions 
starts in an inclusion. Hence, the crack has to exit the inclusion first 
which causes an increase of the effective crack energy. 

Fig. 6.11b and Fig. 6.12b show the scatter of the effective crack energy 
computed for all 100 microstructure realizations in the lower fiber count 
range. For both volume fractions, we observe that the Dirichlet boundary 
conditions result in a much wider range of possible values for the effec- 
tive quantity Yer than the periodic boundary conditions. Furthermore, 
we observe a division of this wide range into wide scatter, about one 
third of the data for 30% filler fraction and on half of the data for a filler 
fraction of 50%. Furthermore, we see an accumulation of the remaining 
data around lower effective values. Additionally, we notice that the 
range of the outliers decreases for increasing fiber count since these 
effects result from initial cracks inside of inclusions. 

To gain additional insight into the influence of the boundary conditions, 
we investigate the median as well as the upper and the lower percentile 
ranges in Fig. 6.11c and Fig. 6.12c. We observe that the range of effective 
crack energies for the two boundary conditions under consideration 
overlap for both volume fractions. Hence, for some microstructures, the 
Dirichlet boundary conditions result in the same effective crack energy as 
the periodic boundary condition for a possibly different microstructure 
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Figure 6.12: Comparison of the boundary conditions for 50% filler content. (Ernesti et al., 
2023) 
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realization. Furthermore, we notice that for both boundary conditions 
the total range and the mid percentile range become smaller for an 
increasing fiber count. The median lines are approaching each other 
as the fiber count increases, however, at a very low rate. In general, 
the range of possible values for the effective crack energy for Dirichlet 
boundary conditions is much wider compared to periodic boundary 
conditions. Furthermore, the median for periodic boundary conditions 
is roughly placed in the center of the data set. In contrast, for Dirichlet 
boundary conditions the median and the mid percentile range are placed 
in the lower quarter of the data sets, reflecting the aforementioned 
outliers. 

Last but not least, we investigate the relative standard deviation of 
the two data sets over the fiber count, see Fig. 6.11d and Fig. 6.12d. 
We observe a decrease of the standard deviation as the fiber count 
increases. Furthermore, we notice that the standard deviation for 
periodic boundary conditions is more than one magnitude lower than 
for Dirichlet boundary conditions. Specifically, for a volume fraction 
of 50%, we notice an increase of the standard deviation for the last 
microstructure sample with 80? fillers. A possible explanation for this 
effect may be found in Fig. 6.12c where we notice that the spread of the 
standard deviation is caused by some computations with lower effective 
crack energy compared to the median/mean value. These lower outliers 
are caused by sections of straight crack paths which are still possible and 
probable for very large microstructures. 

To sum up, we strongly discourage using the Dirichlet boundary condi- 
tions. Rather, periodic boundary conditions should be preferred. 


6.5 Conclusion 


In this chapter we studied the influence of the boundary conditions on 
the effective crack energy of heterogeneous materials. Based on homoge- 
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nization result (Braides et al., 1996; Friedrich et al., 2022; Cagnetti et al., 
2019) for the Francfort-Marigo model of brittle fracture (Francfort and 
Marigo, 1998) in a quasi-static setting and without crack irreversibility, 
we investigated a method for computing the effective crack energy 
using the fast marching method (Sethian, 1996). We validated our 
approach and compared it to the FFT-based methods using periodic 
boundary conditions discussed in chapters 4 and 5. In addition to 
periodic boundary conditions, the fast marching method provides 
additional freedom in the boundary condition choice. With this freedom 
at hand we compared periodic and Dirichlet boundary conditions for a 
continuously reinforced composite with tough inclusions, containing 
filler fractions of 30% and 50%. We noticed in a study with several 
realizations of volume elements of increasing size that the periodic 
boundary conditions result in a much lower spreading of the results 
compared to Dirichlet boundary conditions. This was reflected in the 
standard deviation, which was one magnitude lower for the periodic 
boundary conditions compared to Dirichlet boundary conditions. For 
an increasing size of the computational cell, we noticed that the medians 
approached each other. However, periodic boundary conditions should 
be preferred over Dirichlet boundary conditions due to the much 
lower standard deviation. This lower standard deviation indicates that 
the necessary computational cell for periodic boundary conditions is 
considerably smaller than for Dirichlet boundary conditions. Thus, we 
strongly recommend using periodic boundary conditions. 

Applying periodic boundary conditions in the context of the fast 
marching method relied on an iterative process over one axis of the 
domain, i.e., increasing the complexity of the algorithm on an N x N grid 
from O(N? log N) for Dirichlet boundary conditions to O(N? log N). For 
microstructures of moderate size, i.e., up to N = 2048, the fast marching 
method is still competitive with an FFT-based solver for the minimum 
cut/maximum flow problem. However, for larger structures the higher 
complexity forms a strong argument against using fast marching for 
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periodic boundary conditions. 

Classical fast marching algorithms are only applicable to isotropic 
crack resistances in the plane. To cover anisotropies in the crack 
resistance (Ernesti and Schneider, 2022), anisotropic fast marching 
methods Waheed (2020) may be explored. 

Last but not least, let us mention that it would be desirable to have 
mathematical results at hand which concern the influence of boundary 
conditions on the effective crack energy. Indeed, for elastic solids, 
results (Sab, 1992; Bourgeat and Piatnitski, 2004; Owhadi, 2003) are 
available which provide a list of suitable boundary conditions whose 
influence becomes negligible when going to the infinite-volume limit. 
Previous work by Bouchitte and Suquet (1991; 1994) for limit-load 
problem suggests that Dirichlet boundary conditions may be used, 
whereas Neumann boundary conditions give rise to different results. 
Further research may be necessary to clarify this issue. 
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Chapter 7 


Summary and conclusions 


The two main goals of this thesis were to find suitable characterizers 
for microstructures applicable to a wide class of microstructures, and to 
establish computational tools for multi-scale brittle fracture. 

Well established multi-scale homogenization methods for hardening 
type materials (Matouš et al., 2017) face difficulties when applied to 
softening type materials (Gitman et al., 2007). Therefore, a first attempt 
to simulate phase-field fracture on heterogeneous microstructures, dis- 
cussed in section 2.3, provides only limited applicability for a multi- 
scale approach to fracture. This problem motivated us to pursue a 
different path. 

In this work we focussed on the Francfort-Marigo model of brittle 
fracture (Francfort and Marigo, 1998). The Francfort-Marigo model 
provides a variational approach to fracture by subsequently solving free- 
discontinuity problems fulfilling an irreversibility constraint. Braides 
et al. (1996) provided a periodic homogenization result for a similar 
class of free-discontinuity problems, which is directly applicable to a 
single increment of the Francfort-Marigo model under specific loading 
conditions, namely anti-plane shear, while neglecting the irreversibility 
constraint. Furthermore, extensions lifting these requirements individu- 
ally were established. Giacomini and Ponsiglione (2006) included the 
irreversibility constraint and Cagnetti et al. (2019) provided an extension 
to random, ergodic media, both within the anti-plane shear setting. 
Friedrich et al. (2022) lifted the restriction to anti-plane shear for periodic 


191 


7 Summary and conclusions 


microstructures and the two-dimensional case. 
Motivated by the homogenization result of Braides et al. (1996), which 
gives specific formulas for the homogenized bulk and surface term in 
free-discontinuity problems, the effective crack energy is defined as the 
minimum cut through the microstructure of varying crack resistances. 
To solve this problem, Schneider (2020) proposed an FFT-based solution 
strategy for complex three dimensional microstructures based on the 
minimum cut/maximum flow duality shown by Strang (1983). 
Chapter 3 was devoted to the characterization of digital microstructures. 
Microstructure characterization is concerned with finding simple quan- 
titative expressions to classify microstructures and provides a crucial 
step before finding effective material parameters using homogenization 
approaches. By taking an approach based on Minkowski tensors which 
originate in stochastic geometry, we established the quadratic normal 
tensor as a suitable characterizer for microstructures which, due to local 
symmetries on the microscale, indicates a macroscopic anisotropy in the 
effective material behavior. Our main findings are listed below. 

e We established the quadratic normal tensor as a robust characterizer of 
microstructures since it is translation invariant and does not scale with 
the size of the microstructure. Due to its tensorial form it naturally 
accounts for anisotropy. 

e We provided a robust computational method to compute the quadratic 
normal tensor on large scale gray images of microstructures. We care- 
fully validated our approach and investigated multi-grid convergence. 

e We investigated complex microstructures of industrial interest. First, 
we compared the quadratic normal tensor with the fiber orienta- 
tion tensor of second order for a fiber reinforced composite. We 
furthermore compared its accuracy with a well established method to 
compute fiber orientation tensors, namely a structure tensor approach. 
Finally, we investigated sand-binder structures used in casting appli- 
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cations, demonstrating the wide range of possible shapes for which 
the quadratic normal tensor may be computed. 


In chapter 4 we established a novel discretization and solver strategy for 
the cell formula for computing the effective crack energy. The approach 
of Schneider (2020) used trigonometric collocation (Moulinec and Suquet, 
1994; 1998) and finite elements with reduced integration (Willot, 2015a) 
to discretize the cell formula. Furthermore, an FFT-based primal-dual 
hybrid gradient method (Chambolle and Pock, 2016) was used to solve 
the discretized equations. This initial work had two shortcomings. The 
involved discretizations showed checkerboard artifacts and non-smooth 
solution fields. Furthermore, the performance of the solver prevented 
finding high accuracy solutions. In (Ernesti et al., 2021) we found that 
higher accuracy and smoother solution fields may be established using 
a discretization presented by Couprie et al. (2011) based on a combi- 
natorially consistent grid, namely the CCMF discretization. However, 
only computations using a small number of degrees of freedom were 
possible with this approach as only direct solvers were available. Our 
main findings in chapter 4 were the following. 

e We provided an implementation of the CCMF discretization into a 
framework for FFT-based micromechanics. The implementation relies 
on a doubling of the degrees of freedom. The resulting solution 
fields exhibit far less artifacts than the discretizations previously 
investigated by Schneider (2020). 

To solve the discretized problem we relied on the alternating di- 
rection method of multipliers (ADMM) (Glowinski and Marrocco, 


1975; Gabay and Mercier, 1976). In particular, augmenting the strat- 
egy with an adaptive parameter strategy (Schneider, 2021b) showed 
promising results. 

e Ina study on a continuous-fiber reinforced composite we found that 
for the CCMF discretization the residual correlates with the error 
compared to an accurate solution (Ernesti et al., 2021; Domahidi 
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et al., 2013) which was not the case for other discretizations (Willot, 
2015a). The adaptive strategy showed a strong reduction of the 
iterations required to reach desired accuracies. In particular, we 
found the Barzilai-Borwein adaptive strategy (Xu et al., 2017) to be 
most promising. 

e In subsequent numerical experiments we investigated complex mi- 
crostructures ranging from fiber reinforced composites to structures 
with porous space. 

We provided an extension of the approach presented in chapter 4 in 

chapter 5 to account for locally anisotropic crack resistances. Our main 

findings are listed below. 

e We considered an anisotropic minimum cut problem which includes 
the anisotropy via a tensorial crack resistance. For this case we derived 
an anisotropic maximum flow cell formula to compute the effective 
crack energy. 

e In order to include this anisotropy into the framework presented in 
chapter 4 we established a projection operator, more precisely the 
projection onto an ellipsoid. To solve this projection problem we 
relied on Newton’s method. 

e This novel framework including anisotropy allows the study of ad- 
ditional classes of materials. We investigated polycrystalline brittle 
materials with a cleavage plane in each grain, as well as a carbon fiber 
reinforced composite with a transversely isotropic crack resistance in 
each fiber. 

In chapter 6 we focussed particularly on the case of stochastic homog- 

enization where Cagnetti et al. (2019) provided the theoretical back- 

ground. From a mathematical point of view, the qualitative theory 
of stochastic homogenization and (non-porous) linear elastic materials 

is well established (Papanicolaou and Varadhan, 1981; Kozlov, 1978; 

Sanchez-Palencia, 1980). Whereas the effective properties do not depend 

on the boundary conditions, apparent properties, i.e., those evaluated 
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on cells of finite size, are affected by the choice of boundary conditions. 


In particular, a higher accuracy of the apparent properties is typically 


achieved with periodic boundary conditions. In this chapter we inves- 


tigated the influence of the boundary conditions on the apparent crack 


energy evaluated using the cell formula of Braides et al. (1996). Several 


remarks on our main findings are in order. 


Even before their established homogenization result, Braides and 
Piat (1995) showed that using periodic boundary conditions on a 
single periodic cell results in the same effective surface term for free- 
discontinuity problems than considering the infinite volume limit 
using Dirichlet boundary conditions. 

In order to compare the boundary conditions in the stochastic case 
we established the fast marching method (Sethian, 1996; 1999) as 
a tool to compute the effective crack energy on two dimensional 
microstructures. In particular, using fast marching methods provided 
additional freedom in the choice of the boundary conditions compared 
to the approach presented in chapters 4 and 5. 

After a validation of the fast marching approach we investigated 
the influence of the boundary conditions on various ensembles of 
cells and increasing cell size. Our investigations showed that for 
an increasing cell size the influence of the boundary condition de- 
creases. However, Dirichlet boundary conditions resulted in a large 
scatter for different microstructure realizations compared to periodic 
boundary conditions and in particular a standard deviation of one 
magnitude higher. As a result we strongly recommend using peri- 
odic boundary conditions to compute the effective crack energy on 
heterogeneous microstructures. 


Let us put our multi-scale approach for brittle fracture into perspec- 


tive. The foundation for our approach is laid by a mathematical ho- 


mogenization result for free discontinuity problems. Based on this 


result the effective crack energy is defined via a minimium cut problem. 
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We overcame numerical issues, namely checkerboard artifacts in the 
discretization and slow converging solvers, by introducing a novel 
solver and discretization strategy which allowed for higher accuracy 
solutions than previous work (Schneider, 2020). Anisotropic materials 
may be studied by computing the orthogonal projection on an ellipsoid. 
Investigations of the boundary conditions on 2D structures using fast 
marching methods indicated that using periodic boundary conditions 
results in a decrease of the necessary size of the computational domain 
compared to Dirichlet boundary conditions. 

Nevertheless, the presented approach does have its limitations. First 
of all, the original homogenization result holds only for a restricted 
setup of the Francfort-Marigo model, i.e., neglecting irreversiblity and 
assuming anti-plane shear loading. Subsequent work extending the 
original approach (Friedrich et al., 2022; Giacomini and Ponsiglione, 
2006) only lifted these restictions individually but not all at once. 

A second issue lies in the Francfort-Marigo model itself and its math- 
ematical roots in free-discontinuity problems. The homogenization 
result proves I’-convergence of the functionals upon homogenization. By 
definition, [-convergence implies the convergence of global minimizers 
of the functionals but provides no statement on local minimizers. Fol- 
lowing physical intuition, local minimizers may be more realistic but are 
excluded from the mathematical treatment. Consider for instance a crack 
propagating through a microstructure. As the crack hits an obstacle, it 
may change its direction based on what is currently energetically more 
favorable, but may not follow the globally minimal path. 

A further consequence of the given homogenization result is the de- 
coupling of the effective crack energy and the effective stiffness upon 
homogenization. This may indeed contrast with physical intuition as 
pointed out by Michel and Suquet (2022), who considered a two phase 
laminate material with different elastic properties in the phases but 
equal crack resistance y. Michel and Suquet (2022) pointed out that the 
laminate material may in fact be toughened due to the presence of a 
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phase with a lower stiffness. In this case the effective crack resistance 
may exceed y which is not possible with our approach. 

An additional issue is caused by the presence of two small scales, i.e., the 
characteristic length of the microstructure, and the loading increment on 
the macro-scale, which is considered quasi-static. The ansatz proposed 
here considers a fixed discretization of the load into small increments 
and vanishing size of the microstructure in each step. Therefore, if a 
macroscopic crack propagates by a certain increment as a result of an 
increased load, the size of the microstructure upon homogenization is 
smaller than the crack increment. Hence, the macroscopic crack passes 
the microstructure cell within a single load increment. This point of 
view contrasts with other approaches from the literature on multi-scale 
fracture mechanics. For instance, Hossain et al. (2014) investigated quasi- 
static crack propagation on the microscale using a phase-field fracture 
model and identified the effective crack resistance as the maximum J- 
integral in time. However, their approach is not based on a mathematical 
homogenization result, and thus, a macroscopic model would have to 
be postulated. More to the point, the limit of vanishing time step size 
and the limit of vanishing microstructure size do not seem to commute. 
Hence, further investigations on this may be required. 

Furthermore, investigations whether the quasi-static assumption holds 
may be of interest. In particular, quasi-static refers to a low velocity of an 
applied load which allows to neglect dynamical effects. It also assumes 
that the crack itself travels with a low velocity. If the crack passes the 
microstructure within a single load step its velocity may in fact be too 
large for this assumption. This would require a novel homogenization 
result which accounts for dynamic effects. 

Lastly, let us point out that the Francfort-Marigo model does not distin- 
quish tension and compression since the bulk energy of the Francfort- 
Marigo model is quadratic in the strain field. This contradicts physical 
intuition, see Fig. 2.2 when a compressive load causes the crack phases 
to close which may result in an interpenetration of the material with 
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itself and motivated energy splittings for phase-field fracture mod- 
els (Amor et al., 2009; Miehe et al., 2010b). The model by Amor et al. 
(2009) converges to the Francfort-Marigo model with the additional 
constraint to non-interpenetrating crack faces for vanishing phase-field 
width (Chambolle et al., 2018). For this constraint to non-interpenetrating 
crack faces no homogenization result is yet available. 

To conclude, the presented approach does have its limitations. On the 
one hand, the basis is laid by a proven homogenization result. But 
on the other hand, several assumptions for the homogenization result 
to hold may contradict with physical intuition. To overcome these 
issues, introducing these physical constraints into novel homogenization 
approaches may serve as a first step. From the application point of 
view, further investigations on the non-commuting small scales may 
be conducted. Michel and Suquet (2022) suggested for instance an 
additional time stepping on the microscale. However, caution has 
to be taken in order to stay on the trail provided by a well defined 
homogenization result. This distinguishes our approach from others 
who compute the effective crack resistance (Lebihain et al., 2021; Hossain 
et al., 2014). Furthermore, due to the global minimization involved, 
the established effective crack energy does form a lower bound to 
the true effective crack resistance. It is therefore suitable to model 
the worst case scenario, but may not incorporate certain toughening 
effects (Michel and Suquet, 2022). 

Future work to continue our ansatz may discuss the up-scaling aspect of 
the presented homogenization approach. Using the presented tools one 
may compute the effective crack energy and use well known methods to 
compute the effective stiffness (Milton, 2002) of given microstructures. 
Using these effective quantities, simulations on the component scale 
may be conducted using phase-field methods. These face two main diffi- 
culties in the presence of anisotropy. First, the use of energy splittings 
to account for a tension-compression anisotropy faces difficulties if the 
stiffness is anisotropic. Due to the anisotropy, the eigensystems of local 
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strains and stresses do not align. In this case, an additive decomposition 
of the energy is not possible, in general. Furthermore, the question arises 
how to define tension or compression in the anisotropic setting. An 
initial contribution was proposed by van Dijk et al. (2020) who consider 
the elastic energy as an isotropic tensor function of a quantity Y via 
e:C:e=W:W. For an energy decomposition expressed in Y, the 
well-established splittings of Miehe et al. (2010b) and Amor et al. (2009) 
are applicable. 

Secondly, the phase-field model has to account for an anisotropic crack 
resistance. Several approaches have been proposed to incorporate tenso- 
rial crack resistances, we refer to section 6.1 for an overview. Following 
Focardi (2001) we may chose an approximation of the surface energy of 
the form 


[alt a] dx (7.1) 
a2L! 

within a phase field model for some anisotropic norm p : RT > Rso and 
mean crack resistance y. A common approach (Clayton and Knap, 2014) 
uses y?(Vd) = Vd- MVd for some symmetric, positive definite matrix 
M. One question of particular interest is how to relate the anisotropic 
norm y and the prefactor 7 with the efective crack energy 7°". 

Finally, a strategy for components with locally varying fibrous mi- 
crostructure may be conducted. Here, the approach of fiber orientation 
interpolation proposed by Köbler et al. (2018) may be promising. 
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Minkowski tensors for specific 
shapes 


A.1 Minkowski tensor of a ball 


Consider the ball Bp(0), parameterized by spherical coordinates (r, y, 0), 
with r € [0,R), p € [0,27], and 0 € [0,7]. The transformation to 
Cartesian coordinates reads 


r sin(@) cos(y) 
x(r,y,0) = | rsin(#) sin(y) 
r cos(0) 


The outward-pointing unit normal on 0BR(0) is given by n(r, y, 0) = 
x(1,y, 0) and is thus independent of r. With this parameterization at 
hand, the Minkowski tensor WP’? computes as 


4 7 R? 20 T P 
WIR) = I f n(2,6) 8 n(v, 8) sin(6)dddy 


u AnR? 


9 Id. 
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A.2 Minkowski tensor of a cylinder 


We consider a cylinder K in R3, oriented in z-direction. We parameterize 
it by cylindrical coordinates (r, p,z) with r € [0,R),p € [0,27] and 
z € (0, L). The transformation to Cartesian coordinates reads 


We divide the boundary into three subsets, describing the side, top 
and bottom of the cylinder 0k = OK, U 0K; U OK». The side OK, is 
parameterized by r = R, y € (0,27), z € [0, L], the bottom OK, by 
r € [0, R], y € (0, 27], z = 0 and the top 0K; by r € [0, R], y € (0, 2m], 
z = L. The outward-pointing unit normals for the side, bottom and top 
boundary, respectively, read 


cos(p) 0 0 
n, = | sin(y) |, u =| 0 and m = 0 
0 1 —1 


With this parametrization at hand, we compute the Minkowski tensor 
WP? of K by 


Bu cos? (4) cos(p)sin(p) 0 
wW? (K -3/ a cos(p ‘i sin? (y) 0 | Rdypdz 
0 


2 27 
+- [ / e, ®e,rdrdy 
3 Jo Jo 


2 
= w 8 er + ey Q ey) + = Re, Qez 


L 
= 7 fe- Qe; + — IR (14-0: ae.) 
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A.2 Minkowski tensor of a cylinder 


Dividing WP? (K) by its trace gives the quadratic normal tensor of K: 


QNT (ic) = — @e,+ ao (1 —e, Q e+) : 
If R < L holds, then R/(R + L) is the smallest eigenvalue, which 
indicates an extension in e,-direction. The larger eigenvalue L/(2(R+L)) 
has multiplicity 2, indicating some symmetry within the ey — e,-plane. 
If R > L holds, the smaller eigenvalue has multiplicity 2, indicating a 
disc-like shape within the e, — e,-plane. 
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Performance of additional penalty 
factor choices for ADMM 


In addition to the lower bound and the Barzilai-Borwein strategy for 
choosing the penalty factor p in combination with different damping 
parameters 0, see Section 4.4.3, we investigated two additional choices 
which are popular in the literature. More precisely, we consider resid- 
ual balancing (He et al., 2000) and the averaging strategy proposed 
by Lorenz and Tran-Dinh (2019), which perform admirably for linear 
elastic and inelastic homogenization problems (Schneider, 2021b). The 
resulting residual and error plots are shown in Fig. B.1. For the CCMF- 
discretization and the damping parameter 6 = 0.5, the residual balancing 
strategy led to an unstable behavior. The choice 6 = 0.25 resolves this 
instability. However, this approach does not lead to a highly accurate 
solution. The averaging strategy by Lorenz and Tran-Dinh (2019) shows 
more promising results, reaching a tolerance of 1074 in fewer than 2000 
iterations and ö = 0.25. However, this parameter choice turns out to 
be inferior to the Barzilai-Borwein approach. The relative error (4.36), 
shown in Fig. B.1b correlates with the residual in a similar way as for 
the choices considered in Section 4.4.3. For the rotated staggered grid 
discretization, the Lorenz-Tran-Dinh scaling with ô = 0.25 shows the 
best performance. However, only low accuracy in terms of the relative 
error (4.36) may be reached. 
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Figure B.1: Residual and error measure for CCMF and rotated staggered grid discretiza- 
tions, comparing different solver parameters. (Ernesti and Schneider, 2021) 
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